Yan Rolland

CT Computed Tomography FOV Field Of View FVM Finite Volume Method Gy Gray HCC HepatoCellular Carcinoma LAVA Liver Acquisition with Volume Acceleration LBM Lattice Boltzmann Method MRI Magnetic Resonance Imaging PET Positron Emission Tomography RF Radiofrequency ROI Region Of Interest SIRT Selective Internal Radiation Therapy SNR signal to noise ratio SPECT Single-Photon Emission Computed Tomography TACE Trans Arterial Chemo Embolisation LIST OF

Keywords:

Arrivée à la fin de ce doctorat, je veux remercier toutes les personnes qui ont participé, de façon plus ou moins directe, à cette étape si importante pour moi. Par contre, je m'accorderai ici la liberté de faire des phrases aussi longues que je le voudrai, après tous les efforts de clarté nécessaires pendant la rédaction du manuscrit, je m'en excuse d'avance.

Je ne peux que commencer par ma directrice de thèse, Johanne Bezy-Wendling, qui m'a toujours soutenue et m'a énormément appris aussi bien sur le plan personnel que scientifique. Elle me connaît maintenant tellement bien qu'elle peut répondre à mes questions avant même que je ne les pose. Si ma thèse s'est toujours passée dans la bonne humeur, même pendant les périodes difficiles, c'est surtout grâce à elle et à la bonne ambiance qui a toujours régné au sein de toute l'équipe Metriq.

I would like to thank my co-supervisor, Marek Kretowski, for his insightful and always pertinent advices that improved the quality of my thesis and often helped me finding my way. Thank you very much to Krzysztof Jurczuk too, for all the times he was there to help, I couldn't even count them. Thank you to both of you and to all the people I met in Białystok for the warm welcome during the 7 months I spent in Poland.

Un énorme merci à mon co-encadrant Pierre-Antoine Eliat, qui m'a ouverte au monde de l'imagerie médicale déjà avant le début de ma thèse, lorsqu'il m'a présenté l'IRM 4,7 T du labo et a pu observer le vide dans mes yeux... heureusement les choses ont beaucoup changé depuis et c'est en particulier grâce à lui.

Je tiens à remercier Hervé Saint-Jalmes qui a dépensé tant d'énergie à m'enseigner le fonctionnement de l'IRM, merci pour tout le temps et la patience lors des acquisitions.

comprenais qu'un mot sur 5, c'est grâce à lui si maintenant je me sens prête à injecter des microsphères dans l'artère hépatique.

Ma thèse n'aurait pas été la même sans l'aide fondamentale de Jean-Jacques Bellanger, que je ne remercierai jamais assez d'avoir sacrifié un grand nombre d'après-midis de grandpère épanoui pour venir au labo nous aider et nous donner ses conseils essentiels et décisifs. Ces années n'auraient surtout pas été pareilles sans tous les collègues et amis doctorants qui arrivaient et qui partaient, mais qui venaient toujours au ru à 11h36 et qui avaient toujours quelque chose d'intéressant à raconter. Tout d'abord Hélène, Caroline et Amandine, mes repères omniprésentes respectivement au labo de Beaulieu et à celui de Villejean, qui étaient toujours là et sur qui j'ai toujours su que je pouvais faire confiance pour quoi que ce soit. J'aurais tellement de choses à dire mais ça sortirait largement de ce cadre.

Merci à mes collègues de bureau: Mahmoud qui m'a toujours donné les meilleurs conseils, avec qui j'ai partagé d'innombrables discussions philosophiques et qui me manquera sûrement. Merci à Yvan pour avoir été le plus fiable des fournisseurs de gâteaux et chocolats, alsaciens ou non, et à Nahed pour sa présence rassurante.

Merci à Juan, Juan-David, David, Fred et Momo, qui sont partis avant la fin de ma thèse et à ceux qui ont eu l'honneur de participer au costicafé, costofé, pablosti, appelez le comme vous voulez, Matthieu, Nicolas, Pablo, Vivien, Karim, et tous ceux que j'ai déjà nommé. Merci à Mireia, Nadine, Diego, Daniel, Gustavo, Siouar, Axel, Eugenia, Andres, Sahadatou, Nawras, et à tous les autres qui participent à la bonne ambiance du labo.

Grazie a Roberto e a Elena, con cui finalmente posso dire arzigogolare, mietitrebbiatrice, e soprattutto pota! Chez moi on n'aime pas les remerciements formels mais il y a bien deux ou trois choses très concrètes qui ne seraient sûrement pas arrivées sans les injections de motivation constantes de ma famille: je n'aurais pas fait de lycée littéraire, choix qu'a posteriori je referais mille fois pour une infinité de raisons malheureusement mal comprises. Ils m'ont donné la détermination (et le support, grazie Michele!) qui a fait que je n'ai pas abandonné la licence de mathématiques pendant les premières années où entrer dans la salle de cours était pour moi comme débarquer à Pékin. Je n'aurais jamais eu ni l'envie ni les moyens d'aller faire mon master à Paris, qui m'a enfin amené vers ce doctorat. Merci à Robin, qui m'a laaargement aidée à me lancer dans cette experience rennaise et qui tout le long m'a offert son support moral, scientifique, logistique, culinaire,... et un mot spécial pour mon petit neveu adoré Arturo, que tout le monde au labo doit connaître maintenant, grâce aux innombrables photos et vidéos que je n'arrête pas de montrer comme une vieille mamie.

CONTENTS

3.1

For each of the two populations, this table presents mean and standard deviation of ∆C n a (the hepatic veins enhancement between non-injected phase and the arterial phase), ∆C p (the ratio between the non injected-portal phase enhancement and arterial-portal phases enhancement) and ∆C V C (the difference between the contrast uptake on the vena cava before and after the hepatic veins at the arterial phase) . . . . . . . . . . . . . . . . . . . . . . . Le carcinome hépatocellulaire (CHC) est un problème de santé majeur dans le monde, en particulier dans les pays en voie de développement où apparaissent plus de 80 % des cas dans le monde [START_REF] Giunchi | Epidemiology of Hepatocellular Carcinoma[END_REF]. Dans plus de 90 % des cas, le CHC se développe à partir d'une cirrhose hépatique, elle même générée par obésité, alcoolisme ou par la présence d'une hépatite B ou C. Son incidence est très importante en Asie orientale et en Afrique (plus de 20 sur 100000 cas par an), où ce cancer primitif du foie se développe principalement à partir d'une hépatite B ou C [START_REF] Sherman | Hepatocellular Carcinoma: Epidemiology, Risk Factors, and Screening[END_REF][START_REF] Giunchi | Epidemiology of Hepatocellular Carcinoma[END_REF]. Dans les pays d'Amérique du nord et d'Europe de l'ouest l'incidence est de 2-7 sur 100000 cas pas an, et les facteurs de risque principaux sont l'alcoolisme et l'hépatite C. Dans le monde, le carcinome hépatocellulaire est la cinquième tumeur maligne la plus commune chez l'homme et la huitième chez la femme, avec plus de 700000 nouveaux cas par an [START_REF] Giunchi | Epidemiology of Hepatocellular Carcinoma[END_REF].

The table gives the correlation

Les options de traitement du cancer du foie sont limitées. Seulement dans 15 % des cas le patient peut-être traité de façon chirurgicale (transplantation ou résection). 50 % des patients ont accès à des traitements non-chirurgicaux comme l'ablation percutanée par radiofréquences, la chimiothérapie ou la radiothérapie. Des traitements palliatifs sont administrés dans 35% des cas [START_REF] Andreana | Treatment of hepatocellular carcinoma (HCC) by intra-arterial infusion of radio-emitter compounds: trans-arterial radio-embolisation of HCC[END_REF]. En cas d'embolisation de l'artère hépatique, l'apport sanguin du foie est compensé par la veine porte. En présence d'une thrombose de la veine porte, la radioembolisation se révèle une alternative efficace à la chimioembolisation, étant donné ses effets embolisant mineurs, et pour des résultats comparables pour un cancer à un stade intermédiaire. La radioembolisation intra-artérielle se révèle même plus efficace pour une tumeur à un stade avancé [START_REF] Andreana | Treatment of hepatocellular carcinoma (HCC) by intra-arterial infusion of radio-emitter compounds: trans-arterial radio-embolisation of HCC[END_REF].

Le travail de cette thèse se concentre sur la radiothérapie métabolique interne (RTMI), aussi appelée radiothérapie interne vectorisée, qui consiste en l'injection de millions de microsphères marquées à l'Yttrium 90 dans l'artère hépatique, en amont de la tumeur.

L'Yttrium 90 est un pur émetteur de rayonnement Beta, qui a une énergie maximale de 2.27 Megaelectron-Volt et une moyenne de 0.9367 MeV [START_REF] Murthy | Yttrium-90 microsphere therapy for hepatic malignancy: devices, indications, technical considerations, and potential complications[END_REF].

Deux types de microsphères sont utilisées: les SIR-Spheres ® (SIRTeX Medical Limited, North Sydney, Australia) sont des microsphères en résine, leur diamètre varie entre 20 et 60 µm et leur activité, provenant de l'Yttrium 90 présent sur la surface des sphères, de 50 Becquerel (Bq) par sphère. Les Theraspheres ® (BTG, UK) ont un diamètre allant de 20 à 30 µm. Grâce à leur composition en verre, l'Yttrium est mis à l'intérieur de la sphère pendant la fabrication, ce qui donne une activité finale de 2500 Bq par microsphère.

Par conséquent, le nombre de microsphères à injecter pour une même dose est fortement inférieure dans le cas des Theraspheres ® , ce qui rend le traitement plus rapide et moins douloureux pour le patient, en limitant les effets emboligènes.

Le protocole clinique pour l'administration des Theraspheres ® (cf. Fig. 2.7) prévoit actuellement une première acquisition scanner (tomodensitométrie -TDM), ou de réso-nance magnétique (IRM), avec injection de produit de contraste, dans le but d'identifier l'emplacement, l'ampleur et le type de vascularisation des lésions tumorales.

L'acquisition à différents temps à partir de l'injection du produit de contraste permet d'obtenir des séries d'images à différentes phases, apportant chacune des informations notamment sur la vascularisation des lésions. La tumeur et les artères hépatiques apparaissent en hyper signal pendant la phase artérielle, 20-30 secondes après l'injection. Un lavement de la lésion et une prise de contraste de la veine porte et du parenchyme apparaissent ensuite pendant la phase portale.

Le protocole actuel pour l'administration des Theraspheres ® prévoit, environ une semaine avant le traitement, une première artériographie avec un but diagnostique. Lors de cette opération, le radiologue s'appuie sur les angiographies TDM effectuées précédemment et visualise, par de fréquentes injections de produit de contraste, le réseau d'artères restant à parcourir. Le radiologue identifie alors un endroit d'injection optimal et relâche des macros Une optimisation du protocole d'imagerie a également été considérée, dans le but d'extraire le maximum possible d'informations, tout en essayant de minimiser les effets secondaires, dus notamment aux rayonnements. Ainsi, une séquence de contraste de phase (PC) en IRM a été optimisée dans le but de pouvoir mesurer la vitesse du sang dans des petits vaisseaux comme l'artère hépatique, sujets au mouvement respiratoire. Le but de telles mesures est d'utiliser ces informations dans l'initialisation et la validation d'un modèle numérique de l'écoulement du sang et des microsphères (cf. Chapitre 2).

Dans le Chapitre 3 une méthode non-invasive de détection précoce d'un éventuel dépôt pulmonaire des Theraspheres ® est proposée. Dans les chapitres suivants, le modèle numérique est présenté. Une simulation spécifique patient du tissu et de la vascularisation, sains et tumoraux, a été possible grâce à l'initialisation avec des données patients, préalablement segmentées, d'un modèle de croissance vasculaire précédemment développé [START_REF] Jurczuk | In Silico Modeling of Magnetic Resonance Flow Imaging in Complex Vascular Networks[END_REF], Kretowski et al., 2003b] et adapté au problème de la RTMI. Le réseau artériel, segmenté à partir des images cliniques et simulé par le modèle, a ensuite été utilisé comme domaine de simulation de deux méthodes numériques de mécanique des fluides (en anglais CFD, computational fluid dynamics). La simulation de l'écoulement du sang par la méthode des volumes finis et par la méthode de lattice Boltzmann est décrite dans le Chapitre 5. Les mêmes méthodes sont utilisées comme étape préliminaire pour la simulation du transport des microsphères dans le sang, présentée dans le Chapitre 6 de ce mémoire.

Modalités d'acquisition d'images

L'imagerie clinique constitue le point de départ de ce travail de thèse. Le protocole d'imagerie a été analysé et des propositions d'optimisation ont été faites. Le scanner à rayons X est encore très utilisé en clinique à cause de sa plus grande facilité d'accès.

Cependant, l'IRM a une meilleure sensibilité pour la détection des tumeurs hépatiques [Schneider et al., 2006, Oliva and[START_REF] Oliva | Liver cancer imaging: Role of CT, MRI, US and PET[END_REF] Les principales contraintes dans la définition de la séquence IRM concernent d'abord la limitation temporelle due aux temps cliniques très serrés. Ensuite la taille de l'artère hépatique, qui est proche de la limite de la résolution spatiale de l'IRM (de l'ordre du millimètre), et le déplacement du foie dû aux mouvements du diaphragme pendant la respiration, qui peuvent entrainer un déplacement du foie jusqu'à 3 cm.

Une séquence de contraste de phase exploite la différence de la phase des spins en mouvement, induite par l'application de deux gradient bipolaires (Fig.

General introduction

Hepatocellular carcinoma (HCC) is a major healthy issue worldwide. 700000 new cases appear every year, especially in developing countries, where it originates principally from hepatitis and cirrhosis. In western Europe and North America, obesity and alcoholism contribute as well [START_REF] Giunchi | Epidemiology of Hepatocellular Carcinoma[END_REF] This work was realised thanks to a constant and tight cooperation with clinicians.

Indeed, we are aware of the importance of realising mathematical models that are the closest to every patient's peculiarities, and of the difficulties of clinical applications of models that can easily be too abstract if the issue of practical feasibility is underestimated.

Thanks to this strong partnership we could build the foundations of our work on medical imaging.

Even before the development of computational models, image processing tools are applied to the improvement of the current clinical protocol. In particular, the intent of the new protocol proposition concerns first the extraction of more information specific to the single patient. An optimised phase-contrast MRI sequence is proposed in order to compute blood velocity in small hepatic arteries. This is indeed a key value in the initialisation and tuning of every computational model describing patient-specific blood flow. Secondly, a new non-invasive image-based methodology for early shunt detection is proposed. Excluding patients at an early step of the process increases their chance to benefit of an appropriate therapy, without delay. An extensive description of such issues is addressed in Chapters 2 and 3.

The numerical model we propose can be classified into three separate but complementary issues, respectively presented in Chapters 4, 5 and 6. First, the simulation of a patient-specific hepatic arterial tree is proposed. Patient-specific data are thus extracted from the images available in the current protocol. Liver and tumours segmentation is performed on the CT-scan volumes. An accurate segmentation of hepatic arteries is realised on the 3D angiography (Fig. 2.6a). After a mandatory registration step, the model previously developed in our research team is used in order to produce an entire hepatic arterial tree, statistically equivalent to the patient's one and presenting the tumour typical vasculature in the segmented malignancies [START_REF] Bezy-Wendling | A 3D Dynamic Model of Vascular Trees[END_REF], Kretowski et al., 2003b[START_REF] Mescam | A physiologically based pharmacokinetic model of vascular-extravascular exchanges during liver carcinogenesis: application to MRI contrast agents[END_REF][START_REF] Smith | Physical Audio Signal Processing[END_REF], Jurczuk et al., 2014]. All the computational models currently proposed by the literature of the optimisation of catheter tip position concern simple patient-inspired geometries (down to 17 tree leaves) [START_REF] Basciano | Computer modeling of controlled microsphere release and targeting in a representative hepatic artery system[END_REF], Xu et al., 2016].

Knowing the geometry of the entire arterial tree of the patient appears to be necessary, and angiographic imagery can supply this type of information only at a limited resolution level. Typically, we can expect to image only vessels with a diameter larger than 0.5 mm.

Without a finer simulation of the vascular tree, we could not address the computation of particles distribution on the really existing set of terminal branches of the patient's vasculature. The development of patient-specific artificial 3D liver and vasculature is described in Chapter 4.

In order to achieve the simulation of microspheres transport, blood flow has to be addressed as well. This issue is approached in Chapter 5 from two different points of view:

simulations based on the Finite Volume Method (FVM), developed in the Fluent software (Ansys Inc.) are compared to those performed through the Lattice Boltzmann Method

GENERAL INTRODUCTION

(LBM), using a method developed in [START_REF] Jurczuk | In Silico Modeling of Magnetic Resonance Flow Imaging in Complex Vascular Networks[END_REF]. The former is applied only on the segmented part of the vascular tree, whereas the latter can be computed on the synthetically grown tree, thanks to its easy handling of complex geometries and natural parallelisation [START_REF] Jurczuk | Hierarchical Parallel Approach in Vascular Network Modeling -Hybrid MPI+OpenMP Implementation[END_REF].

The last chapter is dedicated to the actual microspheres injection and transport simulation. The ambitions of this part of the work concern the optimisation of the catheter tip placement, but this model will also help optimising the injection velocity and pressure, the number of flushes, the interval between flushes, etcetera. Microspheres transport is also analysed by means of the same two methods presented in Chapter 5 (FVM and LBM), in order to accomplish an exhaustive comparison. Lattice Boltzmann Method appears to be more efficient when the double influence of fluid on microspheres and vice-versa wants to be studied. This appears to be necessary when the ratio between vessel and particles diameter becomes lower, which is the case in the smallest simulated vessels. An envisaged outline considers the coupling of the two methods (cf. Section 4.4).

Chapter 1

Medical context

Liver anatomy and functionalities

Liver is the second largest organ of the human body, the first one being the skin. Liver belongs to the digestive system and has a crucial role in metabolism. Its functionalities include: blood filtration, metabolism of carbohydrates, proteins, fats and foreign chemicals, and storage of vitamins and iron [Guyton and Hall, 2006]. Liver is situated in the abdominal cavity, just below the diaphragm and its larger part is located on the right side of the human body. Its proximity to the diaphragm, the heart and the lungs makes it a relatively mobile organ. Indeed, following the respiratory cycle, this organ can move up to 3 cm. In the average human adult it weights approximately 1.5 kg. It is divided into left and right lobes, the latter occupying approximately 6 times the volume of the former. These two lobes are divided by the falciform ligament. Two other small lobes exist: the caudate lobe and the quadrate lobe (Fig. 1.1). Between them, the hepatic hilum, also known as porta hepatis, is a slit-like opening through which the vena porta and the hepatic artery enter the organ, as well as the bile ducts and the nerves. As shown in Fig. 1.2, liver is further divided into eight segments, which are fundamental references for physicians, for instance during surgery. Liver connective tissue is an assembly of small hexagonal lobules whose size varies between 1 and 2.5 mm, and where all the exchanges between capillaries occur (Fig. 1.3). At the centre of every lobule the central vein brings blood back to the vena cava, while at each of the six corners of the lateral surface of every hexagon is a region called the portal triad. This is essentially defined by the presence of a hepatic arteriole, a portal venule and a bile duct, and is shared by the three neighbouring lobules.

The vasculature of liver is peculiar, due to its double blood supply. Indeed, portal vein furnishes between 70% and 80% of hepatic blood flow, the remaining 30-20% being provided by the hepatic artery. The portal vein brings about 1050 millilitres per minute, additional 300 millilitres flow through the hepatic artery, for a total amount of about 1350 ml/min, which makes approximately the 27 per cent of the total cardiac output [Guyton and Hall, 2006].

The portal vein arises from the gastrointestinal tract and brings all nutrients that will be processed in liver lobules, while the hepatic artery brings oxygenated blood directly from the aorta. Hepatic artery sprouts from the celiac trunk, an artery bifurcating from the descending aorta before the renal arteries. This represents an important anatomical landmark, having the kidneys a very strong need in blood flow: around the 22 per cent of the total cardiac output. Before the bifurcation leading to the proper hepatic artery, the celiac trunk divides into the left gastric artery and next to the splenic artery. The section after this bifurcation is called the common hepatic artery. From the latter then sprout the gastroduodenal, the supraduodenal and the right gastric arteries. Next, this same vessel is called the proper hepatic artery, and it further divides into left and right hepatic arteries, the latter being usually larger than the former, due to the volume proportion between the two lobes (Fig. 1.4).

Concerning the way blood leaves the hepatic tissues, every central vein situated in every lobule is then collected into the three hepatic veins, which reach the vena cava just after they emerge out of the liver. The hepatic segments are also well distinct by the venous vasculature (Fig. 1.2).

The values of velocity and pressure of blood at the entrance of the proper hepatic artery are known in literature: the mean arterial blood velocity at the right hepatic artery is estimated to be around 30 cm/s in [START_REF] Hubner | Hepatic arterial blood flow velocities: assessment by transcutaneous and intravascular Doppler sonography[END_REF], although many different values are available in the literature, varying between 11.6 and 50 cm/s [START_REF] Dyvorne | Quantification of hepatic blood flow using a highresolution phase-contrast MRI sequence with compressed sensing acceleration[END_REF][START_REF] Tziafalia | Echo-Doppler measurements of portal vein and hepatic artery in asymptomatic patients with hepatitis B virus and healthy adults[END_REF][START_REF] Kito | Doppler sonography of hepatic arterial blood flow velocity after percutaneous transhepatic portal vein embolization[END_REF]. Hepatic arterial pressure is less controversial and is evaluated between 90 and 110 mmHg [START_REF] Ficher | The Liver[END_REF][START_REF] Lautt | Hepatic vasculature: a conceptual review[END_REF]. The diameter of the hepatic artery at liver hilum is also very variable, but it usually varies between 4 and 6 mm [START_REF] Da Silveira | Arterial diameter of the celiac trunk and its branches. Anatomical study[END_REF]. Flow precise characteristics further in the hepatic vessels are [START_REF] Martini | Fundamentals of Anatomy and Physiology[END_REF] very poorly documented and can greatly vary between different individuals. Approximate information relying pressure and velocity to vessels diameter in human systemic circuit exists (Fig. 1.5).

It should be noted that the anatomy described above represents the standard pattern, but many kind of anomalies exist in number, size or position of a structure.

Liver cirrhosis and hepatocellular carcinoma

Alcoholism, hepatitis and overweight are the most frequent sources of hepatic cirrhosis, a disease induced by the replacement of the destroyed hepatic cells by fibrous tissue, which eventually obstructs the portal flow by contracting around the vessels [Guyton and Hall, 2006]. Chronic hepatitis B, hepatitis C and cirrhosis are the principal risk factors for the development of a hepatocellular carcinoma. Nevertheless, their respective influence is not uniform within the different regions of the world which are affected by this disease. The incidence of HCC is the greatest in China and Africa, where the major contributor is the presence of hepatitis B [START_REF] Sherman | Hepatocellular Carcinoma: Epidemiology, Risk Factors, and Screening[END_REF][START_REF] Giunchi | Epidemiology of Hepatocellular Carcinoma[END_REF]. Hepatitis C is the principal risk factor for HCC in Japan, while in Western Europe countries and North America the incidence is considered much lower (2-7 per 100000 cases per year) [START_REF] Giunchi | Epidemiology of Hepatocellular Carcinoma[END_REF], but it has been increasing in recent years [START_REF] Sherman | Hepatocellular Carcinoma: Epidemiology, Risk Factors, and Screening[END_REF]. In such countries the development of an HCC is mainly related to alcoholism and hepatitis C.

Throughout the world, HCC is the leading cause of primary liver malignancy. It is the fifth most common malignant tumour in men, the eighth in women, with more than 700000 new cases worldwide per year [START_REF] Giunchi | Epidemiology of Hepatocellular Carcinoma[END_REF].

Even if the portal vein contributes to 80% of hepatic blood supply, a peculiarity of HCC is that it is mainly vascularised by the hepatic artery. Like most common cancers, it stimulates the generation of more vessels in order to have a direct access to the nutrients. This pathological angiogenic process results in a disordered and dense network which presents a high fractal dimension [START_REF] Ledzewicz | Optimal and suboptimal protocols for a mathematical model for tumor anti-angiogenesis in combination with chemotherapy[END_REF]. This arterial hypervascularisation of HCC is the reason why the arterial supply of blood flow within the tumour is 3 to 7 times greater than in the surrounding noncancerous tissue. Within the tumour, the maximum vascular density and the concentration of cells in the tissue are also higher than in the healthy parenchyma, being the growth ratio of cancerous cells higher than the natural one [START_REF] Folkman | Clinical applications of research on angiogenesis[END_REF][START_REF] Brodsky | Vascular density and VEGF expression in hepatic lesions[END_REF]. According to the American Joint Committee on Cancer [START_REF] Edge | American Joint Committee on Cancer (AJCC) cancer staging manual[END_REF], the classification of different stages of the tumour is based on microvascular invasion and tumour size. Four principal stages exist:

• T1: the tumour is isolated and does not present any microvascular invasion;

• T2: multiple lesions smaller than 5 cm appear, or the tumour presents some microvascular invasion;

• T3a: multiple lesions bigger than 5 cm appear;

• T3b: invasion of a major branch of the portal vein or hepatic vein;

• T4: direct invasion of adjacent organs or perforation of the peritoneum.

Treatments of hepatocellular carcinoma

Treatment options for patients diagnosed with hepatocellular carcinoma are limited. For a tumour in an early stage the possibilities are liver transplantation, resection or percutaneous ablation by radio-frequency or microwave [START_REF] Andreana | Treatment of hepatocellular carcinoma (HCC) by intra-arterial infusion of radio-emitter compounds: trans-arterial radio-embolisation of HCC[END_REF]. Intermediate or advanced stage HCC can be treated with chemotherapy, by means of a drug, like Sorafenib, or the administration of chemo-embolization through the vascular path. Radiotherapy can be administered as adjuvant to other treatments at any stage of an HCC. It is estimated that less than 15% of patients are candidates for surgery, 50% for non-surgical treatments and 35% for best supportive care as immunosuppression or local radiotherapy of metastasis [START_REF] Andreana | Treatment of hepatocellular carcinoma (HCC) by intra-arterial infusion of radio-emitter compounds: trans-arterial radio-embolisation of HCC[END_REF][START_REF] Sposito | Comparative efficacy of sorafenib versus best supportive care in recurrent hepatocellular carcinoma after liver transplantation: A case-control study[END_REF].

The major limitation of radiotherapy is the risk of causing radiation-induced liver disease. For example, external beam radiotherapy requires a radiation dose higher than 70

Gy to destroy a solid tumour, while the healthy liver tolerance dose for this therapy is 35

Gy [START_REF] Grouphealth | Clinical Review Criteria. SIRT (Selective Internal Radiation Therapy) Therasphere and SIR Sphere for Unresectable Hepatocellular Carcinoma[END_REF]. The most limiting factors for chemotherapy are the development of drug resistance and its serious side effects [START_REF] Ledzewicz | Optimal and suboptimal protocols for a mathematical model for tumor anti-angiogenesis in combination with chemotherapy[END_REF]. This is why it is important to develop local and regional treatments.

The typical arterial vascularisation of HCC has led the way to a variety of therapies that have improved during the last decades: intra-arterial injection of chemotherapeutic agents, radio-labeled particles, or non-reactive particles, with the only aim of embolisation, i.e. of preventing the blood from delivering nutrients to the cancerous cells.

Trans Arterial Chemo Embolisation (TACE) is thus considered as the gold standard for treating intermediate stage HCC [START_REF] Raoul | Evolving strategies for the management of intermediate-stage hepatocellular carcinoma: available evidence and expert opinion on the use of transarterial chemoembolization[END_REF]. It consists in the infusion of chemotherapeutic agents via the hepatic arterial blood supply, through a percutaneously placed trans-arterial catheter. This is usually followed by selective arterial embolisation of the tumour vascular supply. TACE is shown to induce a tumour shrinkage ranging from the 10% and the 60% of its original volume [START_REF] Carr | Therapeutic equivalence in survival for hepatic arterial chemoembolization and yttrium 90 microsphere treatments in unresectable hepatocellular carcinoma: a two-cohort study[END_REF]. However, post-embolisation syndrome (fever, nausea and pain) can be severe after TACE, therefore this therapy should not be considered in the presence of Portal Vein Thrombosis (PVT) [START_REF] Andreana | Treatment of hepatocellular carcinoma (HCC) by intra-arterial infusion of radio-emitter compounds: trans-arterial radio-embolisation of HCC[END_REF].

A valid alternative is represented by Trans-Arterial Radio Embolisation (TARE): a way to deliver selectively high radiation dose in the area of interest of the liver. It is a form of interstitial-radiotherapy developed for liver tumours, namely a therapy characterised by the placement of the radioactive drug directly into the region of interest. The most frequently reported agents available for this therapy are the Iodine-131 radio-labelled lipiodol and Yttrium-90 loaded microspheres. The first one is strongly limited by the emission of high-energy gamma radiation to which the healthcare personnel and the patient might be exposed. In this case the patient should be hospitalised for several days in a radiationshielded room [START_REF] Andreana | Treatment of hepatocellular carcinoma (HCC) by intra-arterial infusion of radio-emitter compounds: trans-arterial radio-embolisation of HCC[END_REF]. The latter was developed recently to avoid radioprotection problems and the consequent isolation of patients. Indeed, Yttrium-90 is a pure β rays emitter with a half-life of 64.2 hours. Its maximum energy is 2.27 Megaelectron-Volt and the mean energy 0.9367 MeV [START_REF] Murthy | Yttrium-90 microsphere therapy for hepatic malignancy: devices, indications, technical considerations, and potential complications[END_REF].

It has been shown that TARE techniques have far less embolising characteristics than TACE, thus making the treatment possible even for patients suffering from PVT [START_REF] Andreana | Treatment of hepatocellular carcinoma (HCC) by intra-arterial infusion of radio-emitter compounds: trans-arterial radio-embolisation of HCC[END_REF]. Globally, TARE gives equivalent results than TACE for intermediate stage HCC and it presents fewer side effects. For end-stage HCC, TARE therapies are found to have an even better outcome than TACE [START_REF] Andreana | Treatment of hepatocellular carcinoma (HCC) by intra-arterial infusion of radio-emitter compounds: trans-arterial radio-embolisation of HCC[END_REF].

Selective Internal Radiation Therapy (SIRT)

Although its embolising properties are unclear, Selective Internal Radiation Therapy falls within the TARE framework. SIRT consists in the injection of millions of microspheres containing radioactive Yttrium-90 into the hepatic artery, through a catheter introduced in the femoral artery [BTG IM -Therasphere RoW, 2016]. As shown in Fig. 1.6, microspheres flow directly into the liver tumour via its own vascular network and become permanently lodged into its small blood vessels. Consequently, there is a preferential delivery of microspheres to the tumour capillary bed allowing for higher radiation doses to Microspheres continue to emit radiations during several weeks after injection, but radiation intensity decreases quickly to insignificant levels. The half-life of 90 Y is 64.1 hours, and 94% of the dose is emitted in 11 days while it decays to Zirconium-90 [START_REF] Murthy | Yttrium-90 microsphere therapy for hepatic malignancy: devices, indications, technical considerations, and potential complications[END_REF].

SIRT can be delivered through two different kinds of microspheres: SIR-Spheres ® (SIRTeX Medical Limited, North Sydney, Australia) and Theraspheres ® (BTG, UK). They are both loaded with 90 Y but present many differences: the former ones are made of resin and their diameter ranges between 20 and 60 µm, with a mean diameter of 32.5 µm.

Theraspheres ® have a smaller mean diameter, between 20 and 30 µm, and are made of glass. Their activity is 2500 Becquerel (Bq) per sphere [START_REF] Murthy | Yttrium-90 microsphere therapy for hepatic malignancy: devices, indications, technical considerations, and potential complications[END_REF], 50 times greater than the resin-based spheres. This is due to the different techniques used for their production: for SIR-Spheres ® the Yttrium-90 is situated on the surface of the sphere, whereas in the case of Theraspheres ® it is incorporated into the glass under its nonradioactive form: Yttrium-89. Afterwhile, the microspheres are bombarded by neutrons in a reactor, becoming Yttrium-90 loaded [START_REF] Carr | Hepatic arterial 90Yttrium glass microspheres (Therasphere) for unresectable hepatocellular carcinoma: interim safety and survival data on 65 patients[END_REF]. Moreover, unlike SIR-Spheres ® , TheraSpheres ® do not appear to have embolisation effects. Therefore, this therapy can be used even in case of portal vein thrombosis and it keeps future treatment options open as the patient progresses through the later stages of the disease. This is due to the higher specific activity of glass microspheres, which implies a lower volume of microspheres to dispense the same activity. Therefore, embolic occlusion of the parent artery is not observed arteriographically. Another advantage of glass microspheres concerns the duration of the intervention: a typical glass microspheres injection lasts 2-5 minutes, whereas in the case of resin microspheres it can last up to 35 minutes. Furthermore, due to the high quantity of needed resin microspheres, SIR-Spheres ® injection appears to be more painful for the patient [START_REF] Sangro | Survival after yttrium-90 resin microsphere radioembolization of hepatocellular carcinoma across Barcelona clinic liver cancer stages: a European evaluation[END_REF].

The current protocol for the treatment of HCC includes a first CT-scan to identify and locate the malignancies. Nowadays, the use of MRI is also increasing thanks to its less invasive nature. Both imaging modalities are performed with contrast agents and at different phases of its distribution within the vasculature. These acquisitions allow a proper evaluation of tumour extent, help the radiologist to guide the angiography, and are a baseline reference to assessment of response. The administration of SIRT requires of course a good knowledge of the arterial tree's branches and bifurcations. For this purpose an arteriography is realised through the injection in the liver arteries of a iodine-based contrast media and the consecutive acquisition of a 3D rotational Cone Beam CT-scan (Fig. 2.6a). This acquisition gives a detailed portrait of the arterial tree and it allows the interventional radiologist to choose the injection point in the arteries for future delivery of the treatment.

The choice of this point is currently validated by the injection of a particular kind of human proteins, which are likely to imitate the behaviour of the microspheres in its transport by blood [START_REF] Garin | Clinical impact of 99mTc-MAA SPECT/CT-based dosimetry in the radioembolization of liver malignancies with 90Yloaded microspheres[END_REF]. Such protein is a Technetium 99m Tc albumin macroaggregate ( 99m Tc-MAA). Its distribution is captured by a SPECT-CT acquisition, which also seeks to detect a significant uptake of the aggregate outside the liver. In particular, a pulmonary planar scintigraphy is acquired in order to estimate hepato-pulmonary shunt.

This diagnostic arteriography also permits to determine the dosimetry, in order to calculate the right activity to inject during treatment. Indeed, based on the SPECT-CT acquired one hour after the injection, the nuclear radiologist measures the volumes of left and right lungs, the volume of all the liver tissue where the activity exceeds a certain threshold (the injected liver volume) and the volume of the tumour (Fig. 1.7). Next, based on the activity per pixel, he is able to measure the dose that every tissue has absorbed based on the simplified accepted formula [START_REF] Garin | Clinical impact of 99mTc-MAA SPECT/CT-based dosimetry in the radioembolization of liver malignancies with 90Yloaded microspheres[END_REF]:

D = A(1 -S) 50 W , (1.1)
where D represents the dose delivered to the injected liver (measured in Gray), A the total

The activity to be injected was calculated conventionally for the schedule to achieve an injected liver dose of 120 6 20 Gy (without exceeding 30 Gy for the lungs). Microsphere injection was performed 8-15 d later. The corrected doses were also calculated after treatment administration, taking into account the activity actually injected, IA cor (IA 2 residual activity remaining in tubing).

Response Evaluation

Tumor response was assessed using the criteria of the European Association for the Study of the Liver (8). Triphasic CT scans were obtained 3 mo after treatment and then every 3 mo until progression or death. Response was evaluated according to the European Association for the Study of the Liver as follows: complete response (CR), the absence of any enhancing tissue; partial response (PR), greater than 50% decrease in enhancing tissue; progressive disease (PD), an increase greater than 25% in the size of one or more measurable lesions or the appearance of new lesions; and stable disease, a tumor response between PR and PD.

Statistics

Activities are expressed in GBq, absorbed doses in Gy, and quantitative values as mean 6 SD.

A Wilcoxon distribution-free comparison test was used to compare the doses absorbed by the tumors. The threshold dose (TD threshold ), that is, the tumoral dose that has to be achieved to observe a morphologic response, was determined using receiveroperating-characteristic analysis and the Youden index. To avoid missing a patient who might respond and therefore benefit from treatment, the threshold was defined as having a sensitivity of 100%.

True-positives are defined as objective response (CR 1 PR) obtained when tumoral dose (TD) is higher than or equal to the TD threshold . True-negatives are defined by a TD less than TD threshold , without an objective response at 3 mo. False-negatives are cases in which TD is less than TD threshold and shows an objective response at 3 mo. False-positives are cases in which TD is at least as great as TD threshold and shows no objective response at 3 mo.

Univariate analysis and a Wilcoxon distribution-free comparison test were used to identify parameters associated with tumor response, progression-free survival (PFS), and overall survival (OS). Data selected after univariate analysis were then subjected to multivariate analysis using a logistic regression test. PFS is the time between treatment and progression (CT scan) or death or second cancer. OS is the time between treatment and the last follow-up visit or death. The survival rates were estimated using the Kaplan-Meier method, and their comparison was based on a log-rank test.

We used SAS software (SAS Institute Inc.) for the statistical analyses, with a significance threshold set at a P value less than 0.05.

RESULTS

Fifty injections of 90 Y-loaded glass microspheres were performed for 36 patients. Three patients did not receive any therapy (2 because of high hepatopulmonary shunting and 1 for failure of the angiography procedure). Fourteen patients received 2 treatments (bilateral lesions, 3; incomplete injection of the initial dose, 2; poor response after the first procedure, 4; and recurrence after treatment, 5). Fiftyeight lesions were assessed.

Table 2 shows the data concerning IA (planned, IA plan , and corrected, IA cor ), biodistribution, and dosimetry. The mean value of the injected healthy liver dose was 79.9 6 24.5 Gy (minimum,19.6 Gy;maximum,194 Gy), repre- MAA SPECT/CT DOSIMETRY IN LIVER • Garin et al. 257 by on October 3, 2016. For personal use only. jnm.snmjournals.org Downloaded from

The activity to be injected was calculated conventionally for the schedule to achieve an injected liver dose of 120 6 20 Gy (without exceeding 30 Gy for the lungs). Microsphere injection was performed 8-15 d later. The corrected doses were also calculated after treatment administration, taking into account the activity actually injected, IA cor (IA 2 residual activity remaining in tubing).

Response Evaluation

Tumor response was assessed using the criteria of the European Association for the Study of the Liver (8). Triphasic CT scans were obtained 3 mo after treatment and then every 3 mo until progression or death. Response was evaluated according to the European Association for the Study of the Liver as follows: complete response (CR), the absence of any enhancing tissue; partial response (PR), greater than 50% decrease in enhancing tissue; progressive disease (PD), an increase greater than 25% in the size of one or more measurable lesions or the appearance of new lesions; and stable disease, a tumor response between PR and PD.

Statistics

Activities are expressed in GBq, absorbed doses in Gy, and quantitative values as mean 6 SD.

A Wilcoxon distribution-free comparison test was used to compare the doses absorbed by the tumors. The threshold dose (TD threshold ), that is, the tumoral dose that has to be achieved to observe a morphologic response, was determined using receiveroperating-characteristic analysis and the Youden index. To avoid missing a patient who might respond and therefore benefit from treatment, the threshold was defined as having a sensitivity of 100%.

True-positives are defined as objective response (CR 1 PR) obtained when tumoral dose (TD) is higher than or equal to the TD threshold . True-negatives are defined by a TD less than TD threshold , without an objective response at 3 mo. False-negatives are cases in which TD is less than TD threshold and shows an objective response at 3 mo. False-positives are cases in which TD is at least as great as TD threshold and shows no objective response at 3 mo.

Univariate analysis and a Wilcoxon distribution-free comparison test were used to identify parameters associated with tumor response, progression-free survival (PFS), and overall survival (OS). Data selected after univariate analysis were then subjected to multivariate analysis using a logistic regression test. PFS is the time between treatment and progression (CT scan) or death or second cancer. OS is the time between treatment and the last follow-up visit or death. The survival rates were estimated using the Kaplan-Meier method, and their comparison was based on a log-rank test.

We used SAS software (SAS Institute Inc.) for the statistical analyses, with a significance threshold set at a P value less than 0.05.

RESULTS

Fifty injections of 90 Y-loaded glass microspheres were performed for 36 patients. Three patients did not receive any therapy (2 because of high hepatopulmonary shunting and 1 for failure of the angiography procedure). Fourteen patients received 2 treatments (bilateral lesions, 3; incomplete injection of the initial dose, 2; poor response after the first procedure, 4; and recurrence after treatment, 5). Fiftyeight lesions were assessed.

Table 2 shows the data concerning IA (planned, IA plan , and corrected, IA cor ), biodistribution, and dosimetry. The mean value of the injected healthy liver dose was 79.9 6 24.5 Gy (minimum,19.6 Gy;maximum,194 Gy), repre- [START_REF] Garin | Dosimetry Based on 99mTc-Macroaggregated Albumin SPECT/CT Accurately Predicts Tumor Response and Survival in Hepatocellular Carcinoma Patients Treated with 90Y-Loaded Glass Microspheres: Preliminary Results[END_REF] activity (in GigaBecquerels) to be injected, S the percentage of the hepatopulmonary shunt and W the mass of the treated liver (kg). The formula is based on the principle that 1 Gigabecquerel (GBq) of 90 Y delivers 50 Gy to a 1 kg lesion. Such computation allows the physician to measure the pulmonary shunt and the dose to be injected. The former must not exceed 30 Gy for a single administration, or 50 Gy for multiple administrations. The latter is based on the criteria that the dose delivered to the injected liver should be 120±20 Gy, although [START_REF] Garin | Boosted selective internai radiation therapy with 90 Y-loaded glass microspheres ( B-SIRT ) for hepatocellular carcinoma patients : a new personalized promising concept[END_REF] have proposed a method of dose personalisation where the dose delivered to the tumour is "boosted" up to 205 Gy, resulting in an amelioration of time to progression and of overall survival.

Nevertheless, healthy liver tissue tolerable dose in SIRT is still discussed in the literature. The higher tolerable dose in SIRT compared to external beam radiotherapy is probably related to the non uniformity of the microspheres microscopic distribution, allowing for tissue regeneration [START_REF] Högberg | Increased absorbed liver dose in Selective Internal Radiation Therapy (SIRT) correlates with increased sphere-cluster frequency and absorbed dose inhomogeneity[END_REF]. Indeed, [START_REF] Kennedy | Pathologic response and microdosimetry of (90)Y microspheres in man: review of four explanted whole livers[END_REF] and [START_REF] Högberg | Increased absorbed liver dose in Selective Internal Radiation Therapy (SIRT) correlates with increased sphere-cluster frequency and absorbed dose inhomogeneity[END_REF] among others showed the presence of microspheres clusters in arterioles and small arteries. This means that while some regions of the tissues receive a very high radiations, other are almost spared and, due to the high regeneration capacity of liver cells, they help healthy tissue regrowth [START_REF] Cremonesi | Radioembolization of hepatic lesions from a radiobiology and dosimetric perspective[END_REF]. The case study described by [START_REF] Högberg | Increased absorbed liver dose in Selective Internal Radiation Therapy (SIRT) correlates with increased sphere-cluster frequency and absorbed dose inhomogeneity[END_REF] concerns however a patient treated with resin microspheres, and no precise data are furnished for the clustering formation of glass microspheres. [START_REF] Kennedy | Pathologic response and microdosimetry of (90)Y microspheres in man: review of four explanted whole livers[END_REF] displayed the presence of clusters also in the case of glass microspheres, still noticing the higher number of microspheres for resin clusters, due to the higher number of needed microspheres, for the same dose, compared to the glass ones.

One to two weeks after this diagnostic arteriography, the therapeutic arteriography delivers the microspheres loaded with 90 Y and another SPECT/CT is acquired to observe microspheres distribution. Efficacy is finally assessed by morphologic and functional imaging (either MRI or CT-scan) performed every 2 to 3 months, with objective response commonly occurring 4 to 6 months after treatment. Several parameters are currently used to estimate the efficacy of anti-cancer treatments. Some are simply based on tumour size (WHO or RECIST criteria [START_REF] Therasse | New guidelines to evaluate the response to treatment in solid Tumors[END_REF]), others on vascular modifications, which are shown to occur sooner [START_REF] Bruix | Clinical management of hepatocellular carcinoma. Conclusions of the Barcelona-2000 EASL conference. European Association for the Study of the Liver[END_REF]. Some criteria reflecting vascular changes have later been introduced (mRECIST, [START_REF] Lencioni | Modified RECIST (mRECIST) Assessment for Hepatocellular Carcinoma[END_REF]). Functional imaging techniques reflecting vascularity and perfusion modifications, as well as tumour cells changes due to treatment, have also been proposed but not standardised yet [START_REF] Kim | Assessment of the treatment response of HCC[END_REF]. Examples include diffusion weighted MRI [START_REF] Atassi | Multimodality Imaging Following 90Y Radioembolization: A Comprehensive Review and Pictorial Essay[END_REF] and MRI with new contrast agents like Gd-EOB-DTPA.

Current limitations to SIRT

As outlined before in this chapter, liver vasculature is very complex and only an invasive arteriography can currently help the physician during the planning and the administration of SIRT. Indeed, physicians have to tackle multiple difficulties:

1. the choice of the injection point;

2. the assessment of potential extra-hepatic uptake; 3. the evaluation of the exact activity to be injected. Currently, all the three points are handled through the previous administration of 99m Tc-MAA, leading the patient to undergo two invasive arteriographies within one or two weeks.

The choice for the placement of the catheter tip in the hepatic artery is made by the radiologist without any kind of technological help, making it difficult to reproduce and subject to human error. Its only validation is the 2D and, if needed, 3D CT-scan acquired during such diagnostic intervention, and the scintigraphy acquired a posteriori.

The injection of 99m Tc-MAA is also currently fundamental to determine if the patient can be accepted into the SIRT protocol or not. Indeed, the hepatopulmonary shunt is measured based on the scintigraphy post 99m Tc-MAA injection. The risk of extra-hepatic uptake of radioactive microspheres, namely to the lungs or the digestive tract is the principal contraindication to SIRT. The former might appear in the case of an arteriovenous hepatic shunt: if the blood finds a way to reach the hepatic veins directly from the arteries, without passing through the capillaries, the microspheres will not be captured in the small vessels and will reach the lungs. A cumulated radiation dose to the lungs exceeding 50 Gy is an absolute contraindication to SIRT [START_REF] Andreana | Treatment of hepatocellular carcinoma (HCC) by intra-arterial infusion of radio-emitter compounds: trans-arterial radio-embolisation of HCC[END_REF]. The latter might occur due to a reflux through the gastroduodenal artery, which arises from the common hepatic artery as shown in Fig. 1.4.

The activity to be injected is also based on the dosimetry measured on the distribution of the 99m Tc-MAA, but actually it can happen that the physician has to interrupt the 90 Y injection earlier than planned, due to some resistance in the vessel. It is important to notice that every prevision concerning glass microspheres based on MAA is controversial, since their behaviour can be even very different, due to their different size and nature. It is currently impossible to have a precise prevision of the distribution of microspheres when transported by blood. Even the scintigraphy acquired after the injection can only provide an approximation of their distribution map, that can be erroneously interpreted due to the acquisition system. For instance, 90 Y is much less captured by the bremsstrahlung SPECT/CT scan used in clinical routine than by other available scans like 90 Y PET/CT [START_REF] Elschot | Quantitative Comparison of PET and Bremsstrahlung SPECT for Imaging the In Vivo Yttrium-90 Microsphere Distribution after Liver Radioembolization[END_REF][START_REF] Kao | Post-radioembolization yttrium-90 PET/CT -part 1: diagnostic reporting[END_REF]. Another limitation to SIRT is its high cost compared to chemoembolisation or to radiofrequency ablation [START_REF] Ray | Interventional radiologic treatment of hepatocellular carcinoma-A cost analysis from the payer perspective[END_REF].

Globally, selective internal radiation therapy is more efficient than other HCC therapies [START_REF] Garin | Boosted selective internai radiation therapy with 90 Y-loaded glass microspheres ( B-SIRT ) for hepatocellular carcinoma patients : a new personalized promising concept[END_REF], but it still needs to be optimised in terms of time and cost, as well as personalisation. A numerical model of liver tissue and vasculature, as well as blood and microspheres flow, would help the physician in the choice of the best injection site and dose amount, with the aim of tumour targeting optimisation.

State of the art about SIRT modelling

The numerical estimation of the distribution of microspheres inside the liver is crucial for the administration of SIRT. A numerical model of the treatment would open the way to the possibility of numerically test different scenarios, analyse different results and ideally automatically identify the optimal conditions for tumour targeting. If the numerical model is accurate enough, the distribution of the microspheres in the liver could be predicted, and the 99m Tc-MAA injection could be avoided, making the protocol much less invasive. The amount of the dose to be injected is also matter of discussion and needs to be optimised. The literature proposes a wide range of mathematical methods contributing to the simulation of SIRT treatment. As shown in Fig. 1.8, numerous steps related to very different fields are needed to achieve a complete simulation integrating every aspect of the optimisation cycle. A state of the art concerning each of the subjects treated in this work will be given at the beginning of the chapters dedicated to the each different kind of applications. In the first chapter, only the models devoted in particular to the simulation of SIRT are described.

The following works will be ordered following the steps proposed in Fig. 1.8, even if a strict classification is not possible due to their frequent combinations.

A simplified simulation of vascular tree and drug distribution and absorption is proposed by [START_REF] Schwen | Spatio-Temporal Simulation of First Pass Drug Perfusion in the Liver[END_REF]. Their recently published model concerns micro and macro vascularisation, as well as drug perfusion. The model considers firstly the simulation of vasculature and tissue of the mouse liver, considering only two vascular trees, the one supplying blood, the other one draining it, and the parenchyma, represented at the length scale of lobules (∼ 2 mm). The algorithm is initialised with two vascular trees segmented from in vivo micro-CT imaging, skeletonised and converted to a binary graph. The tissue being divided in a certain number of lobules, the algorithm aims at developing the initial tree in order to connect to every lobule, based on the criteria of the minimum total blood volume, while keeping a constant blood flow at every lobule. Compounds perfusion in the simulated tree is performed solving the advection equation on the simulated blood vessels.

Last, physiologically-based pharmacokinetic modelling [START_REF] Pang | Advanced Pharmacokinetic Models Based on Organ Clearance, Circulatory, and Fractal Concepts[END_REF] is used for the simulation of the metabolisation of compounds. Limitations of this work include a simple representation of vascular network and of a constant blood velocity. Drug distribution is not simulated explicitly, but it is approximated through the solution of simplified advection equation. [START_REF] White | Building a 3D Virtual Liver: Methods for Simulating Blood Flow and Hepatic Clearance on 3D Structures[END_REF] propose a similar method for the simulation of representative dog hepatic vasculature and tissue. A finer simulation of liver tissue is proposed, but a simplified vasculature is considered, presenting up to 50 bifurcations. The vascular model is developed starting from a 2D initial tree, then transformed into a pseudo 3D model in order to avoid vessels overlap between the two trees. Similarly to [START_REF] Schwen | Spatio-Temporal Simulation of First Pass Drug Perfusion in the Liver[END_REF], the vascular growth method is inspired by the work of [START_REF] Schreiner | Computer optimization of vascular trees[END_REF], who proposed the constrained constructive optimisation (CCO). The algorithm is also based on an initial vascular tree whose growth is governed by the pressure values given to randomly positioned leaves, iteratively connected to the tree while minimising the intravascular total volume. Blood flow and drug metabolism are developed through the finite-difference solution of convection-diffusion-reaction equations as the Darcy's law, solved on both vasculature and tissues, considered as porous media. This approach propose an even more simplified simulation of a planar vascular network than [START_REF] Schwen | Spatio-Temporal Simulation of First Pass Drug Perfusion in the Liver[END_REF], presenting non realistic values of blood pressure and velocity.

A similar vascular growth model as the one described by [START_REF] Schwen | Spatio-Temporal Simulation of First Pass Drug Perfusion in the Liver[END_REF] is used by Walrand et al. [2014a,b] for the simulation of SIRT, the optimisation criteria for tree growth being the minimisation of the total vessels length. A Monte Carlo particle transport method is proposed in order to describe the non uniformity of the dose absorption in the tissues, at the level of micro-vasculature and in lobules. Microspheres distribution from larger arteries until portal triads is simulated by taking into account a constant symmetric or asymmetric probability at every bifurcation. Mathematical simulation of microscopic dose distribution is performed based on the Russell's dose deposition kernel. In this model, precise simulation of microspheres transport is not performed, assuming a fixed distribution probability.

The literature also proposes models aiming more specifically at the simulation of SIRT treatment. The more detailed one is proposed by [START_REF] Kleinstreuer | A new methodology for targeting drug-aerosols in the human respiratory system[END_REF], who originally focused on the simulation of the distribution of drug-aerosols inhaled in the respiratory system. Such model was further enhanced and applied to hepatic tumour targeting in SIRT. [START_REF] Basciano | Computer modeling of controlled microsphere release and targeting in a representative hepatic artery system[END_REF], [START_REF] Kennedy | Computer Modeling of Yttrium-90-Microsphere Transport in the Hepatic Arterial Tree to Improve Clinical Outcomes[END_REF] proposed a computational model able to simulate the transport of resin or glass 90 Y loaded microsphere in a representative hepatic artery. Computational Fluid Dynamics (CFD) equations are solved on a domain designed in SolidWorks, which imitates the typical shape of the first 3 bifurcations of the common hepatic artery. In this domain, a CFD software developed by Ansys Inc. is used to solve, through the finite volumes method, the governing equations where blood variable viscosity is modelled by the simplified Quemada model [START_REF] Buchanan | Rheological effects on pulsatile hemodynamics in a stenosed tube[END_REF]. Particles transport is treated by an Eulerian-Lagrangian approach. Fluid values at every node of the mesh are needed for the solution of the Newton's 2nd Law, which governs particles trajectories. The forces considered in such phenomenon are the frictional and drag force, fluid pressure and gravity. Inter particles collisions are neglected due to the dilute nature of the suspension, while particle-wall interactions are defined through the definition of the normal and tangential restitution coefficient (cf. section 6.2.1). The influence of the catheter tip radial position on microspheres distribution was studied. In particular, a particle release map is drawn on the cross section of artery where the catheter tip is placed. The backtracking of particles trajectories allows to show the influence of the injection radial position on the exit branch of the particle. Such method was validated by the use of a physical phantom around 4 times bigger than its anatomic equivalent in [START_REF] Richards | Experimental microsphere targeting in a representative hepatic artery system[END_REF].

The use of a specified micro-device to position the catheter in the vessel section is later proposed by the same research team. The influence of such catheter presence on the flow is studied in [START_REF] Kleinstreuer | A new catheter for tumor targeting with radioactive microspheres in representative hepatic artery systems. Part I: impact of catheter presence on local blood flow and microsphere delivery[END_REF], and a patient-inspired geometry is introduced in [START_REF] Childress | A new catheter for tumortargeting with radioactive microspheres in representative hepatic artery systems. Part II: solid tumor-targeting in a patient-inspired hepatic artery system[END_REF]. The influence of the injection moment, during the cardiac cycle, on particles distribution was studied in Childress and Kleinstreuer [2014a], comparing a full transient simulation of the cardiac cycle to multiple steady simulation considering the different phases. Due to the strong computational time difference and the lack of significant precision improvement, a time-average steady simulation is suggested. Similar conclusions are drawn in Childress and Kleinstreuer [2014b], when a simulation taking into account vessel wall flexibility is compared to a simulation launched on a rigid-diastolic geometry and on a rigid-time averaged geometry. Recently, a more developed patient inspired planar artery geometry with 17 outlets was used [START_REF] Xu | Solid Tumor Embolotherapy in Hepatic Arteries with an Anti-reflux Catheter System[END_REF]. In this work, as well as in [START_REF] Umbarkar | Computationally efficient fluid-particle dynamics simulations of arterial systems[END_REF], more attention has been devoted to the effects of pressure drop at the outputs of the simulated domain on blood flow, due to catheter deployment and embolisation induced by particles presence. Nevertheless, glass microspheres are shown to have a negligible embolising capacity, so this last model does not represent the behaviour of Theraspheres ® but only of resin spheres. The method proposed by this research team is limited by the compulsory use of the proposed micro-catheter in order to select the exact injection point in the vessel section. The concrete application of this method from in vitro to clinical routine appears to be problematic in all centres proposing SIRT treatment and can appear as a restraint to a personalised treatment.

Moreover, no methodology aiming to extract patient-specific vessels geometries in clinical routine is illustrated, and the proposed geometries are limited in number of bifurcations and affinity to reality. Aramburu et al. [2016a] propose a CFD blood flow simulation in the patient-specific hepatic artery, reconstructed from clinical images, using a lumped parameter model for the definition of output conditions. More details will be presented in Section 5.1.

Last, two works describing only the radiations distribution in the tissue will be presented. [START_REF] Gulec | Hepatic Structural Dosimetry in 90Y Microsphere Treatment: A Monte Carlo Modeling Approach Based on Lobular Microanatomy[END_REF] proposed a model of lobular micro-anatomy based on simple geometrical shapes. In this work a periodic microspheres distribution in the tissue is considered (with a differentiation between healthy tissue and tumour), and the deposition of β energy from microspheres in the tissue is simulated through a Monte Carlo model. Different microstructural dosimetry model are compared, showing the much higher radiations concentration at the portal tracts compared to surrounding parenchyma and central veins.

The project of [START_REF] Petitguillaume | OEDIPE , a software for personalized Monte Carlo dosimetry and treatment planning optimization in nuclear medicine: absorbed dose and biologically effective dose considerations[END_REF] focuses on the following three steps: first, the simulation of a patient specific numerical 3D voxelised phantom, based on the patient's CT-scan or MRI. The phantom is realised by means of the in-house developed software OEDIPE [START_REF] Petitguillaume | OEDIPE , a software for personalized Monte Carlo dosimetry and treatment planning optimization in nuclear medicine: absorbed dose and biologically effective dose considerations[END_REF]. Secondly, a mathematical estimation of the curve of the activity distribution over time, based on multiple scintigraphies acquired at different times. Thirdly, a Monte Carlo patient specific simulation of the energy deposition in every organ or tissue, as a function of the delivered activity. Energy deposition is indeed something that can not be measured experimentally, and depends on many factors like tissue density and heterogeneity. The commonly used formula (1.1) is indeed only one of the possible simplifications.

The last two works do not consider neither a hepatic vasculature model, nor the simulation of the injection of radionuclides in the arterial tree.

Problem statement

As illustrated in Fig. 1, numerous problems have to be solved in order to apply the SIRT therapy efficiently. This requires to take into account the peculiarities of a given patient, starting from the extraction of its anatomic characteristics from the medical images, until the computation of the microspheres distribution in the thinnest vessels, perfusing its tumour and healthy tissues.

The goal of this thesis is to contribute to the implementation of the workflow illustrated in Fig. 1.8, taking into account the patient's characterisation through the image processing step, and focusing on the modelling steps, whose different options are illustrated in Fig. 1.9. Such model must be based on a geometrical quantification of the principal part of the hepatic arterial tree, and on tumour regions identification through segmentation/ skeletonisation / registration processing of images, deriving from different acquisition modalities (presented in Chapter 2), the description of the image processing methods and of the resulting workflow being postponed to Chapter 4. The interest of blood velocity measures (cf. Fig. 1.9), led us to phase contrast (PC) MRI acquisitions, whose the principle and experimental results are presented in Chapter 2.

Modelling steps, an overview of which is given in Fig. 1.9, are investigated in Chapters 5 and 6.

9"

8"

2" "

1"

1'" 3" CFD0"

4" " CFD1" «"One" way"»" "

3'" CFD0"

4'" CFD1" «"One" way"»" 5'" CFD2" «"Two" ways"»"

Par;al"problem:"Par;cles"transport" in"the"pa;ent"main"tree"(extracted" from"images)"

Complete"problem:"Par;cles" transport"in"the"pa;ent"main"tree" extended"by"small"simulated"vessels" 6" Boundary"condi;ons" 1D#simplified#approach:#" simula;on"of"a"complete" vascular""tree" (structure"+"flow)" Inputs:"" 1)"pa;ent's"main"tree" (extracted"from"images)" +"tumour"and"liver"shapes" 2)"Input"and"output"pressures" and"flow"(measures"+"literature)"" Outputs:"simplified"geometry" (tubes)""of"the"complete"tree"+" simplified"computa;on"of"flow" and"pressure"by"Poiseuille's"law." 0D" Velocity" measurements" CFD":"FVM"V"LBM" 7" Flow and transport can be simulated through CFD only on the principal arterial tree, segmented from patient images data (on the left of the red line in Fig. 1.9) and closed by suitable boundary conditions, or on an extended tree (on the right), obtained by extending the segmented one, through an image data guided simulation.

Box 7 in Fig. 1.9, untitled "simplified 1D simulation" corresponds to an extension of the principal tree in a "data-fidelity" manner, that is proposed in Chapter 4. The constraints to be taken into account in this context were (cf. Fig. 4.1):

• the patient's principal tree, obtained from cone-beam angiography,

• the patient's liver and tumour(s) contours, obtained from CT-scan images.

The ambition was then, considering that blood flow and pressure at the inlet and at the arterioles level, in tumour and healthy tissues, have to be given as complementary inputs, to obtain as outputs of this simplified simulation:

• the diameter of every vessel branch,

• the blood flow in every vessel branch,

• the pressure drop along every vessel branch.

In order to limit the computational times, a chosen constraint was to derive the results from the Poiseuille's law.

In Fig. 1.9, the arrow going from Box 7 to Box 6 indicates that flow and pressure information, derived from the simplified simulation, can be used as BC in a CFD simulation in the simple principal (segmented) tree. Distal pressure conditions are difficult to obtain in the literature and at the same time are known to have a crucial influence on the simulations [START_REF] Grinberg | Modeling blood flow circulation in intracranial arterial networks: A comparative 3D/1D simulation study[END_REF].

Box 2 indicates that in Chapters 5 and 6, two types of CFD are considered: the classic FVM (Finite Volume Method) simulation and the LBM (Lattice Boltzmann method) simulation. Such kind of simulation has been considered interesting in particular when the microspheres size is non negligible, relatively to the vessel size. In particular, the LBM is classically used for the simulation of flow in complex and tortuous domains, with complex borders as porous media, and therefore appeared to be well suitable for the simulation of blood flow in a complex hepatic arterial tree. An important advantage of the LBM is its natural parallelisation, that has already been verified in the research team [START_REF] Jurczuk | GPU-based computational modeling of magnetic resonance imaging of vascular structures[END_REF][START_REF] Jurczuk | In Silico Modeling of Magnetic Resonance Flow Imaging in Complex Vascular Networks[END_REF]. Abas et al. [2016a] realised a comparison of numerical performance between LBM and FVM in a comparable context: blood flow is simulated in a simple vessel geometry, with a larger diameter than hepatic artery. It is shown that the advantage of LBM against FVM decreases as long as the number of cores increases. Nevertheless, the Lattice Boltzmann method still appears to be more efficient. Concerning the blood flow and microspheres transport simulations, the symbols CFD0, CFD1 and CFD2 (appearing in Boxes 3, 4, 3', 4', 5' ) correspond, respectively, to:

• simple flow simulation (CFD0), which is acceptable if it is supposed that the particles follow accurately the streamlines of the fluid,

• computation of particles trajectory taking into account Newton's second law, without altering the flow lines (one way computation, CFD1)

• the two ways approach (CFD2), which takes into account the fact that the fluid streamlines are more or less modified when the particles are transported.

Box 8 indicates that velocity measures, acquired in this thesis from MRI when available in some branches of the vascular tree, can be taken into account for the model validation and for BC.

To complete the description of the Fig. 1.9, Box 9 concerns the 0D simulation approach, which reduces the CFD problem to the computation of voltage drops and intensities in an electrical discrete components network. Such approach can be used in certain cases to estimate unknown BC [Aramburu et al., 2016a].

Chapter 2

Image acquisition modalities

All the principal steps of this thesis are based on the analysis of medical images and the consequent extraction of patient's data. The numerical simulation of blood flow and microspheres transport, which is our final target, is performed on a computational domain obtained from the images realised in the context of the hepatocellular carcinoma imaging protocol. The initialisation of such kind of hemodynamic models is also a major issue in the literature, due to the serious lack of studies providing blood pressure and velocity values in small vessels like hepatic arterial network. Imaging techniques allow the detection and extraction of a wide range of information while being for most cases totally non-invasive, or minimally invasive. Our attention was then naturally led to an investigation of the current clinical protocol and to an improvement of it in order to optimise and enlarge the possible set of data obtained from the patient's images.

This chapter provides a description of the current imaging protocol applied for patients suffering from HCC and candidate for a Theraspheres ® treatment, next to a brief description of each of the involved imaging techniques. Next, section 2.2 will report about the work devoted to the improvement of the current protocol, through the introduction of some MRI sequences or the revision of the already adopted ones.

Imaging techniques of the current clinical protocol

A general description of each of the involved imaging techniques is given before the actual clinical imaging protocol, applied in the Theraspheres ® framework, will be detailed.

Ultrasonography

Ultrasound (US) imaging is a totally non-invasive modality based on the reflection of sound waves sent towards the considered organ. A probe is used to send the sound wave and receive its echoes reflected by the tissues. Different densities of the tissues produce different strength of the echo they send back, and thus different contrasts on the image.

The time elapsed between the moment the wave is sent and the moment its echo is detected indicates the distance between the tissue and the transducer. This modality is not frequently used in liver cancer protocol, but it can help in characterisation of fibrosis and identification of tumour stage. Nevertheless, it is considerably less accurate than MRI or CT-scan and can not detect lesions smaller than 1 cm. Doppler effect is used to evaluate the movement of the considered structures. A practical application is the study of blood flow. Doppler ultrasounds allow to determine if blood is flowing towards or away from the transducer and it can measure its relative velocity.

Concerning the analysis of cancer, portal flow study can help in quantifying changes in hepatic blood supply. Nevertheless, the accuracy provided by such modality is lower than, for example, X-ray angiography realised with the injection of contrast product. Moreover, despite its low cost and easy accessibility, the quality of the information provided by ultrasonography is strongly dependent on the operator and on the patient's characteristics.

In our study, we realised Doppler US in order to validate blood velocity measurements performed with MRI (cf. section 2.2.2.4).

X-ray Computed Tomography

Commonly called CT-scan, X-ray Computed Tomography is an anatomical imaging technique based on the absorption of an X-rays beam by the tissues. A source of X-rays is placed over the field of interest and a multi-detector is positioned on the opposite side.

The source rotates around the patient, and the acquired signal is processed in order to reconstruct the structure of the organ from its projection on the detectors. The contrast between different tissues is produced by their different absorption coefficients. This technique can be used in order to produce 2D or 3D images. In the first case, a fan-shaped beam of X-rays is projected on the multi-detector which receives one-dimensional information (Fig. 2.1, left). In this way, a stack of 2D images is acquired, in order to reconstruct an image volume. An intrinsic 3D technique is the Cone-Beam CT-scan. As the name [START_REF] Scarfe | Clinical Applications of Cone-Beam Computed Tomography in Dental Practice[END_REF] suggests, a cone-shaped X-rays beam is directed toward a two dimensional array of detectors (Fig. 2.1, right). Fewer rotations of the source are needed in this case, reducing the acquisition time.

CT-scan offers a good image resolution, usually under the millimetre, while allowing entire liver imaging during the time of one breath hold, thus reducing movement artefacts.

Intravenous contrast agent (mostly iodine-based) is often used in order to enhance the image contrast. Furthermore, contrast agent distribution in liver vessels and tissue over time gives fundamental information for the characterisation of the cancerous nodules. As shown in Fig. 2.2a, before the injection tumour and vessels have lower signal intensities than hepatic parenchyma. Approximately 25-30 seconds after the injection, contrast agent reaches the liver through the hepatic arteries, and is still not spread in all the tissues. This phase is called arterial phase and is characterised by the enhancement of the hepatic artery and the tissues mainly vascularised by it, like HCC and metastasis (Fig. 2.2b). The portal-venous phase occurs between 60 and 80 seconds after the injection. In this phase, the injected contrast agent that first went through the gastrointestinal tract, reaches the liver from the path of the portal vein. Maximal enhancement of hepatic parenchyma is seen in this phase (Fig. 2.2c). Concerning the tumour, this is the phase where it is easier to evaluate if it presents an internal necrotic tissue. Indeed, all the hypovascular regions, characterised by a low enhancement, are revealed in this phase.

CT-scan imaging presents a lower cost than MRI and comparable performances in terms of spatial and temporal resolution. Nevertheless, it offers a more limited range of clinical exams and image contrasts. CT-scan is also limited by the need to expose the patient to ionising radiations, and contrast media are usually more nephrotoxic and likely to provoke allergies than the ones used in MRI [START_REF] Namasivayam | Adverse reactions to intravenous iodinated contrast media: A primer for radiologists[END_REF].

Nuclear imaging modalities

This branch of medical imaging is also called functional, physiological or metabolic imaging, due to its faculty to provide important data related to the physiological function of the organ investigated, like metabolic activity of tumour cells. The peculiarity of this modality is that the radiations are emitted from radionuclides previously injected into the patient. Commonly used radionuclides emitting gamma rays for diagnostic nuclear imaging are Technetium-99m, Iodine-123 and Thallium-201. Concerning radionuclides emitting positrons, only Fluorine-18 ( 18 F) is currently used in clinical routine for diagnostic purposes [START_REF] Petitguillaume | OEDIPE , a software for personalized Monte Carlo dosimetry and treatment planning optimization in nuclear medicine: absorbed dose and biologically effective dose considerations[END_REF]. A promising method, used so far for research only, considers the use of Carbone-11 prior to the injection of 18 F [START_REF] Liu | Radiation Dosimetry of Whole-Body Dual Tracer 18F-FDG and 11C-Acetate PET / CT for Hepatocellular Carcinoma[END_REF].

Radioactive isotopes are usually coupled to a carrier (or vector), likely to deposit on a particular biological target. Such administered couple is called a radiopharmaceutical or

radiotracer. An example or radiopharmaceutical used in the Theraspheres ® framework is the 99m Tc coupled to albumin macro-aggregates. The principle of 90 Y loaded microspheres is similar even if they have a therapeutic target and not a diagnostic one.

Two nuclear medicine techniques commonly used in oncology are the Single-Photon Emission Computed Tomography (SPECT) and the Positron Emission Tomography (PET).

In both cases, a gamma camera (also called scintillation camera) is used in order to detect the radiations emitted by the radiotracers from within the patient's body. Image from [START_REF] Soret | Principe et technique de la tomographie d'émission de positons couplée à la tomodensitométrie[END_REF] In the case of PET, gamma rays are pro- should be used [START_REF] Elschot | Quantitative Comparison of PET and Bremsstrahlung SPECT for Imaging the In Vivo Yttrium-90 Microsphere Distribution after Liver Radioembolization[END_REF].

SPECT and PET modalities need to be coupled to standard X-rays CT, due to their poor spatial resolution and to the lack of anatomical contrast. Such hybrid systems called SPECT-CT and PET-CT provide the visualisation of superimposed anatomical and functional information, thus allowing a localisation of the regions with a high concentration of radiopharmaceuticals that could not be possible with nuclear imaging alone.

Angiography

Angiography is an imaging technique dedicated to the depiction of blood vessels. In most cases a contrast agent is injected into the bloodstream for a better image contrast. Angiography can be performed using projectional X-ray radiography, CT or MRI. Cone-beam CT can be used for the realisation of a three dimensional angiography. An example of conebeam CT and radiography realised before the 99m Tc-MAA injection is given in Fig. 2.6a.

MRI is increasingly used, due to its non-invasive nature. Indeed, very detailed images of the vascular network can be performed even avoiding the injection of contrast product.

Magnetic Resonance Imaging

Physical principles of magnetic resonance imaging (MRI) will be described in section 2.2.1.

The use of MRI continues increasing since its first adoption in clinical routine in 1984.

This imaging technique offers a clinical examination comparable to X-rays CT in terms of sensitivity and specificity [START_REF] Oliva | Liver cancer imaging: Role of CT, MRI, US and PET[END_REF]. In particular, MRI has a great sensitivity to image contrasts naturally induced by the different characteristics of tissues. MR images can present many kinds of different contrasts depending on the acquisition parameters and without contrast agents, as will be described in section 2.2.1.

Along with its non-invasive nature, advantages of MRI include better detection and characterisation of tumours than CT. Different branches of magnetic resonance include elastography and spectroscopy: the former reveals mechanical properties of soft tissues, useful for cancer or fibrosis diagnosis [START_REF] Mariappan | Magnetic Resonance Elastography: a Review[END_REF], the latter provides analysis of metabolic changes related for example to cirrhosis or HCC [START_REF] Kuo | In Vivo Proton Magnetic Resonance Spectroscopy of Large Focal Hepatic Lesions and Metabolite Change of Hepatocellular Carcinoma before and after Transcatheter Arterial Chemoembolization Using 3.0-T MR Scanner[END_REF].

Another advantage of MRI is its more flexible character. Whereas CT images can be acquired only on the axial plane, in MRI any of the three principal planes (axial, coronal or sagittal) can be chosen, as well as any oblique one.

Contrast media used in MRI have shown to be much safer than those used in CT-scan and are based on a different principle. Used contrast agents, mostly based on gadolinium, enhance the contrast of the studied tissue by modifying its relaxation time (cf. sec. 2.2.1).

The different possibilities of image contrasts and the injection of contrast media make MRI a particularly suitable technique for the study of liver pathologies [START_REF] Schneider | MRI of the Liver. Imaging Techniques, Contrast Enhancement, Differential Diagnosis[END_REF]. In particular, dynamic contrast-enhanced MRI allows a better characterisation of liver tumours than CT, through the analysis of its modified vascular supply (cf. sec. 2.1.2).

MRI is not limited by the need to minimise the acquisition duration like CT (due to X-ray exposure) and thus allows, for example, multiple acquisitions during the arterial phase of the contrast agent distribution in the liver. In particular, 7 different images are acquired for the follow up of contrast agent in the HCC imaging protocol. As shown in Fig. 2.7, after the first image acquired before the injection, three different acquisitions are realised every 9 seconds during the arterial phase of the contrast distribution. As for CT, another image is acquired during the portal phase when tumour "wash-out" occurs . An equilibrium phase is reached between 3 and 5 minutes after the injection, contrast agent filters into the extracellular space of normal tissue, which continues to be more enhanced than the tumour, where the distribution is mostly intravascular. Two acquisitions are then

The current imaging protocol

A general overview of the current imaging protocol applied for patients candidates for a Theraspheres ® treatment is given in Fig. 2.8. If a patient is believed to be affected by a hepatocellular carcinoma, a CT or MRI examination is first performed. The choice of MRI would generally be preferred due to its non-invasive character and to its better sensitivity and specificity in HCC detection and characterisation [START_REF] Oliva | Liver cancer imaging: Role of CT, MRI, US and PET[END_REF]], but CT is still more frequent due to its lower cost and easier access. This first examination allows the physician to mark the number, the location, the extent and stage of the cancerous nodules.

The arterial vasculature, supplying the tumour, is also analysed thanks to the injection of contrast agent. Due to the importance of MRI in this work, that will be clarified below in this chapter, only the MRI protocol will be described in detail. 5. Between equilibrium and late phases a T2-weighted fat suppression image is acquired to detect and characterise even lesions that may be confused in surrounding fat.

If the tumour stage is too advanced for surgery and if the patient is considered for SIRT, a diagnostic angiography (cf. Fig. 2.8) is performed in order to evaluate the possibility to actually realise the SIRT treatment. It consists in a sort of practical simulation of the treatment, or pre-treatment. During this operation, the interventional radiologist places the catheter into the hepatic artery while being assisted by X-ray radiography of the vessels.

In some cases, the realisation of a 3D cone-beam CT angiography appears necessary for a better overview of the vascular network. The aim of this treatment simulation is to identify an optimal injection position and to inject a radiotracer that is likely to behave similarly to Theraspheres ® . Technetium-99m albumin macro-aggregate are injected and their distribution is captured approximately one hour later by the realisation of a SPECT-CT. This examination provides quantitative data on the distribution of the radiopharmaceutical in the tumour, healthy tissue and even, if any, in surrounding non targeted organs. Moreover, the quantification of 99m Tc-MAA in liver and tumour is used for the planning of 90 Y dosimetry, with the aim of limiting radiations to healthy parenchyma and lungs. Vessels that are shown to bring radiations to surrounding organs may also be coiled.

Between ten days and one week later, the same examination is performed, this time being a therapeutic angiography. The interventional radiologist places the catheter exactly at the same point as for the injection of 99m Tc-MAA and the nuclear physician injects the 90 Y loaded microspheres. A bremsstrahlung SPECT-CT is realised one hour later in order to verify the microspheres deposition.

After the therapy administration, the tumour progression is monitored by regular acquisitions of CT or MRI every 2 or 3 months.

Proposed improvements of the current clinical protocol

The protocol currently used in clinical routine for a patient suffering from HCC still presents scope for improving. X-rays CT is commonly used, exposing the patient to the radiations of contrast product and X-rays themselves. Much information that would help the physician in the characterisation of the tumour, in diagnosis and surgery planning are still inaccessible. An MRI sequence aiming to quantify hepatic arterial blood velocity is presented after a brief, non-exhaustive description of the physics fundamentals of MRI [START_REF] Bernstein | Handbook of MRI pulse sequences[END_REF][START_REF] Kastler | Comprendre l'IRM : manuel d'auto-apprentissage[END_REF]. The detection of hepatopulmonary shunt is possible thanks to the previous realisation of an invasive angiography, where the patient is exposed to X-rays and a catheter is guided along its arteries twice in a short period of time. Replacing this first angiography does not seem to be feasible for the moment, due to the large amount of information it offers for treatment planning. Nevertheless, in Chapter 3 we discuss an alternative way for shunt detection through a modified MRI sequence.

We decided to focus on MRI due to the numerous advantages aforementioned and to the strong background in MRI investigation of the research team. We had the possibility to realise acquisition tests on the clinical 1.5 T MR system (GE Healthcare, Little Chalfont, Buckinghamshire, United Kingdom) of the Eugène Marquis Centre for the fight against Cancer (Centre Eugène Marquis, CEM) of Rennes, France. Images were acquired on test objects filled with a saline solution, healthy volunteers and patients from CEM.

Spatial localisation in MRI and image reconstruction

The principal fundamentals of magnetic resonance are given in Appendix A. Spatial encoding, in particular the one related to the spins phase, needs a more detailed description in view of the sequence we propose and is given below.

In order to create a 2D or 3D image from the signal received by coils, spatial localisation is needed in order to distinguish the origin of the signal in the body. In order to do so, gradients coils are used to apply gradients on the main magnetic field in the three directions, the idea being to slightly modify the character of the signal depending on its position, exploiting the Larmor equation. When the FID signal (Fig. A.2) is detected by the MR system, it fills in, generally line after line, a matrix of raw data called k-space (introduced by [START_REF] Ljunggren | A simple graphical representation of fourier-based imaging methods[END_REF]), to whom an inverse Fourier transform is applied in order to obtain the final image.

1. Slice selection, as well as frequency and phase encoding presented below, can be achieved along any direction. The case of an image acquired on the axial plane is presented. Image slice is selected by the application of a linear magnetic field gradient G z along the z axis. In this way the resonance frequency of the spins becomes proportional to its z coordinate. It is thus sufficient to choose a RF pulse whose frequency bandwidth includes the frequencies of the protons in the slice of interest in order to ignore all the remaining tissues. The application of other gradients in the same way on other directions would give rise to the selection of an oblique slice.

2. The localisation along the x axis is achieved through the application of a frequency encoding gradient G x during the appearance of the echo, namely when the signal is being received by the radiofrequency coils. Since this gradient is applied when the signal is read, it is also called readout gradient. The frequency encoding gradient 2.2. PROPOSED IMPROVEMENTS OF THE CURRENT CLINICAL PROTOCOL 73 is applied in the x direction and it allows to identify the x coordinate of a rotating spin thanks to its frequency during the reading phase. The frequency sampling, realised while the signal is read by the system, is equal to the receiver bandwidth.

The number of sampling steps during the receiver reading time corresponds to the number of columns of the k-space.

3. Phase encoding gradients are applied for the localisation along the last direction:

between the RF pulse and the signal measurement, a gradient along the y axis (G y ) is applied during a short period of time. As a result, after the application of this gradient (also called phase encoding gradient), spins will restart precessing at the same frequency, but they will be out of phase. The different phase will then be the value used to identify the y coordinates of the signal received. Such sequence is applied a number of times corresponding to the number of rows of the k-space, (and thus, generally, of the pixels in the y direction), using slightly different gradients. Every time the elementary cycle is repeated, one new row is filled. This is why the number of pixels in the phase encoding direction strongly influences the acquisition time. A drawback of this method is the phase wrap-around artefacts that can occur due to the periodic character of the phase of spinning protons. Phase wrap-around effects occur if the chosen field of view do not cover the patient in the phase encoded direction. Protons of regions exceeding the field of view have a phase shift φ bigger than the maximum (360°) and will consequently appear folded on the region of interest, superimposed to protons whose phase shift equals φ -360°. Two examples of image artefact due to phase wrapping can be seen in Fig. 2.9, the image field of view does not cover the whole body in the phase encoding direction (here anterior to posterior), thus the excluded posterior part is superimposed on the left of the image.

The acquisition time of a standard MR image depends on the repetition time TR, the number of rows of the k-space and the number of repetitions that may be realised for a single image, for a better signal to noise ratio. MRI allows a large flexibility in the choice of the encoding directions. Axial, coronal or sagittal images can be acquired by a simple exchange of the directions of the encoding gradients. Any oblique plane can be chosen by a combination of the different gradients. the difficulty to perform precise measures on small vessels like hepatic artery moving with diaphragm, as well as the limits in time due to clinical routine.

The phase contrast sequence

The MRI phase contrast (PC) sequence belongs to a category of sequences whose target is the imaging of moving elements. In particular, it is a kind of magnetic resonance angiography (MRA). The PC sequence takes advantage of the difference in the phases of the spins induced by the application of velocity encoding gradients. Two acquisitions are made and then subtracted in order to enhance the signal originating from the moving protons,

PROPOSED IMPROVEMENTS OF THE CURRENT CLINICAL PROTOCOL 75

and to remove the one coming from the static tissues or by unwanted phase differences.

A bipolar gradient, presenting two lobes equal in length and intensity but with opposite sign, noted by Gd + /Gd -, is applied twice in an inverse order. Even if arterial blood flow is pulsatile, PC imaging assumes that blood velocity is constant. Under such hypothesis, it can be assumed that the spins dephasing after the application of a gradient is proportional to their velocity.

During the application of the first bipolar gradient (Gd + /Gd -), the dephasing of the protons of the stationary tissues is first accelerated and then canceled as in a standard gradient echo sequence. The dephasing of moving spins can not be canceled: a proton that is dephased of an angle φ 1 during Gd + , will have moved when Gd -is applied, and it will not be subjected to a -φ 1 rephasing, but to the angle associated to the gradient applied to its new position (for example, φ 3 in Fig. 2.10). The phase shift φ of such moving proton is thus proportional to its velocity and is more precisely described by the equation φ = VγTS,where V represents the velocity of the proton along the direction of the velocity encoding gradient (Gd + /Gd -), γ is the gyromagnetic ratio of the spin, T is Picture taken from [START_REF] Kastler | Comprendre l'IRM : manuel d'auto-apprentissage[END_REF] the time elapsed between the centres of the two lobes of the bipolar gradient, and S is the area under the gradient curve like those presented in Fig. 2.11. In the case described in Fig. 2.10, φ is equal to φ 1 + φ 3 .

A second image is then acquired with the same bipolar gradient, but in the opposite order: Gd -is applied before Gd + . Static tissues behave in the same way, and no dephasing results at the end of this acquisition. The dephasing gathered by the moving spins is the same as during the first bipolar gradient, but with opposite sign: φ = -VγTS.

Two images are acquired thanks to the two acquisitions, presenting an image contrast proportional to the spins phases. The diagram of a typical sequence for the acquisition of one of the two images is shown in Fig. 2.11. In order to enhance the signal relative to flowing blood, the two images are finally subtracted pixel by pixel. The signal deriving from the static tissues is identical in the two images and will be canceled off. The signal deriving from flowing blood is equal but with opposite sign and, as a result of the subtraction, will be doubled.

The maximum dephasing of precessing spins is 180°for spins moving in the positive direction (referring to the first lobe of the first bipolar gradient) and -180°for protons moving in the opposite direction. If the phase shift exceeds this value, protons whose phase difference is 180°+ α could not be distinguished from the protons with a 180°α phase difference. The velocity corresponding to the protons shifted of an angle of 180°i s thus the maximum velocity that can be encoded in the image. This value is called velocity encoding (Venc) and is a parameter that can be chosen by the operator for every acquisition, depending on blood velocity in the studied vessel.

Phase contrast sequences can be used to acquire 2D or 3D images and in any of the three common directions, depending on the direction of the bipolar gradient. As shown in Fig. 2.11, the bipolar gradient associated to one of the three directions (z, x and y)

is added to the gradient used for the spatial localisation in the relative direction (slice selection gradient, frequency encoding gradient and phase encoding gradient). Coupling of the three gradients allows velocity encoding in any arbitrary direction.

Generally, based on the k-spaces of the two acquired images, another image is produced, generally for anatomic orientation. It is called the magnitude image since its contrast is proportional to the magnitude of the spins magnetisation vectors. In most cases it is computed as an average of the magnitude images obtained from the k-spaces of the two acquisitions [START_REF] Bernstein | Handbook of MRI pulse sequences[END_REF].

Limitations of the contrast phase sequence

A frequent clinical application of the phase contrast sequence concerns brain imaging.

The reason lies in the absence of movement related to cardiac or respiratory cycles and to the consequent high precision achieved in such quantitative imaging. Indeed, artefacts related to unwanted protons movement are one of the principal drawbacks of PC sequences.

Concerning liver imaging, ideally it is acquired in a period of time short enough for the patient to be able to hold his breath, typically not more than 20 seconds. If the acquisition is too long (for example in the case of 3D imaging), triggering can be employed, but this can make the acquisition time much longer. If the acquisition is realised during normal respiration, movement artefacts occur and image quality is strongly reduced.

Generally speaking, the phase encoding direction is the most sensitive to movement artefacts. The duration of the phase encoding gradient application is indeed longer than the signal reading, during the frequency encoding gradient application. Moreover, every line of the k-space is filled in after the application of slightly different phase encoding gradient, and this is also likely to induce movement artefacts.

Velocity encoding A difficulty in the correct acquisition of PC images concerns the choice of the velocity encoding value. The challenge is to avoid aliasing without losing too much signal intensity from vessels of slower blood, while encoding the targeted velocities.

Indeed, if the field of view of the acquisition includes protons moving at a higher velocity than Venc, aliasing will occur. This phenomenon is an analogy of the phase repetition: if a proton is moving according to a certain velocity v = Venc + v α , its contrast on the final image could not be distinguished from protons moving at a velocity Venc -v α . An example taken from the images acquired during our study is proposed in Fig. 2.12. Generation of phantoms Phase shifting related to proton movements is exploited for the generation of MR angiography. Nevertheless, if the sequence parameters, mainly TE and TR, are not correctly tuned on the blood vessel that is to be imaged, the uncontrolled phase shift can be a source of artefacts. If the acquisition is tuned on small vessels, repetitions of "phantoms" of vessels with a high blood flow appear along the phase direction. The repetition distance of such phantoms is proportional to the repetition time TR. Indeed, a longer TR means that more time elapses between two elementary cycles, and the protons that are excited during the first cycle have time to go further from the region of interest.

A compromise must then be found in order to obtain an appropriate image quality in a reasonable time. Fig. 2.13 shows an example of aorta and vena cava repetition along the phase encoding direction.

Ground truth Validation of in vivo flow quantification based on PC MRI sequences represents probably the most important challenge. Ground truth is difficult to obtain, especially in hepatic artery. [START_REF] Lotz | Cardiovascular Flow Measurement with Phase-Contrast MR Imaging: Basic Facts and Implementation[END_REF] assess the velocity measurement errors to approximately 10%, while Doppler ultrasonography has been shown to overestimate peak blood velocity of a 25% factor [START_REF] Hoskins | Accuracy of maximum velocity estimates made using Doppler ultrasound systems[END_REF]. Average velocity values at the entrance of the hepatic artery exist in the literature, but no accordance is established. For example, [START_REF] Hubner | Hepatic arterial blood flow velocities: assessment by transcutaneous and intravascular Doppler sonography[END_REF] propose a study based on the comparison between transcutaneous and intravascular Doppler sonography, assessing mean blood velocity in the right hepatic artery [START_REF] Tziafalia | Echo-Doppler measurements of portal vein and hepatic artery in asymptomatic patients with hepatitis B virus and healthy adults[END_REF][START_REF] Hubner | Hepatic arterial blood flow velocities: assessment by transcutaneous and intravascular Doppler sonography[END_REF][START_REF] Kito | Doppler sonography of hepatic arterial blood flow velocity after percutaneous transhepatic portal vein embolization[END_REF]. In those studies blood flow is measured in either the common, left or right hepatic artery, and the average velocity varies between 11 and 50 cm/s.

Choice of other acquisition parameters Thanks to the flexibility in the choice of the direction of the velocity encoding gradients, PC MRI sequences can acquire 2D or 3D

images. Measurement precision in 2D acquisitions is generally improved if the acquisition plane is selected in order to be perpendicular to blood flow. This can lead to acquire images on oblique planes that can be difficult to locate based on standard anatomical references.

A non-negligible issue also concerns the difficulty for the operator to select the correct plane during clinical routine.

Spatial resolution choice has an important influence on image quality. If the chosen voxel size is too small, the number of protons inside a single voxel is little, and the signal deriving from such voxel will consequently be low, leading to a poor signal to noise ratio (SNR). On the other side, if the voxel is too large, different tissues with different characteristics will contribute to the signal deriving from the same voxel, and could not be distinguished. These kinds of drawbacks are called partial volume artefacts.

The challenge concerning hepatic artery imaging is particularly related to spatial resolution. Its small diameter (∼ 5 mm at liver hilum) and its relatively low blood flow would naturally lead to the choice of a small pixel size. Nevertheless, as illustrated in the following section, the SNR can be strongly affected.

Proposed improvements

Acquisitions were performed at the 1.5 T MR system (Optima MR450w, GE Medical Systems) of the Eugène Marquis Centre (CEM), Rennes, France. Three acquisitions were first realised on a test object, then six on healthy volunteers and test object placed under the back of the volunteer, in order to validate the measures. The optimised sequences were finally introduced in the clinical MRI protocol and acquired on 8 HCC patients.

The operator defined parameters that have the most important influence on image quality are velocity encoding value (Venc), voxel size, echo time (TE), repetition time (TR), flip angle (θ), and the receiver bandwidth. The direction to be encoded by the phase or by the frequency gradient can also be modified and has an important role in the development of the artefacts related to phase wrapping.

Acquisitions on test object

Test object consists in a plastic tube connected to the contrast agent injector and, at the end, to a reservoir. The syringes normally used for contrast agent were filled with saline solution and used to inject the solution in the tube at a fixed flow rate. The tube presents three sections characterised by different internal diameters: 1 mm, 2 mm and 3 mm. The injection was launched at flow rates of 0.2, 0.3, 0.4 and 0.5 ml/s. The sequence used in this first part of the study was available in others clinical protocols and not optimised yet, indeed these first tests aimed at validating the velocity measurements performed on the phase image. The sequence parameters are given in Table 2.1. Depending on the MR systems constructor, the conversion between the grey level of the phase image and the actual velocity can be indirect. Some of the constructors introduce an additional process in order to make anatomy more clear on the phase images. In the case of GE Medical Systems, phase images can be magnitude weighted and this information is stored in the DICOM tag called VasFlags. If magnitude weighting is applied to the image, the bit 2 of the vasFlags tag is set.

In this case, the conversion between the image grey level and the velocity in mm/s is given by: encoding value (Venc) is given in mm/s by the dicom tag number (0019,10cc), whereas the velocity encode scale is given by the tag number (0019, 10e2). The vasFlag corresponds to the tag number (0043,1032) and the second bit is set (26=11010), meaning that magnitude weighting is applied. Concerning the 2D phase contrast sequence that we realised, the velocity encode scale is equal to 8.6 and the vasFlag is also set to 26, meaning that magnitude weighting is applied as well. In our experience, the Vscale value was equal to 0.018 for all the 3D PC images, and to 0.054 for all the 2D images, regardless of the chosen Velocity encoding value.

V = Phase pixel
A PC acquisition can measure blood velocity in any of the three directions. A different image is acquired for every direction, and another image is automatically produced by the computation of the norm, pixel by pixel, of the three velocity images. Such image automatically provides the absolute velocity in every direction. All the images reconstructed by the MR system were analysed: the phase image, the magnitude image, and the three images providing the velocity values on the three directions.

3D phase contrast optimisation

Our target being the development of an MRI sequence able to quantify blood velocity as far as possible in the hepatic arterial vascularisation, we first concentrated on the optimisation of a 3D sequence in order to have, ideally, all the required information in one only image.

Velocity encoding values equal to 10, 30, 50 and 100 cm/s were tested. Artefacts related to aliasing were observed in the vena cava, aorta and beginning of the celiac trunk if a Venc ≤ 50 cm/s was chosen.

Slice thicknesses ranging from 1.6 to 8 mm were tested. A very good SNR was observed in the case of a thicker slice (Fig. 2.13), but this induced difficulties in the actual imaging of hepatic artery. Unless the vessel is perpendicular to the acquisition plane, if the slice thickness exceeds the half or the third of the diameter of the vessel to be studied, partial volume effects take place and quantification is compromised. Identifying a vessel by following the bigger one it sprouts from is also problematic if the shift between a slice and the following one is too high and the studied vessel is visible only on one or two slices.

The choice of the flip angle θ and of the repetition time TR are strongly related.

Velocity in each one of the three directions is evaluated in a 3D PC acquisition. If there is one direction that is parallel to the vessel, the same protons can belong to the same slice for a certain time. In this case, if the flip angle is too large and the repetition time is not long enough, the protons that are excited by the θ pulse do not have enough time to recover their net magnetisation, and signal saturation occurs. This is why we defined a small flip angle and a repetition time equal to 10.7 seconds for the 3D images. Echo time was defined as the minimum available on the machine in order to reduce artefacts related to breathing movements. Working on a clinical machine indeed do not allow a complete freedom in the definition of parameters.

Consequences of the modification of the receiver bandwidth are outlined in Fig. 2.15.

As already mentioned, when the signal is detected, a read out gradient is applied, its target being to slightly modify the frequencies of the protons depending on their position along the frequency direction. In this way, the acquired signal includes different frequencies that are centred around a given value depending on the MR system (64 MHz for 1.5 T systems). The width of the frequency domain spectrum depends on the signal sampling rate, on the gradient duration and amplitude. In order to enhance the signal to noise ratio, it is possible to modify the frequency bandwidth of the received signal. A restriction of the bandwidth around the frequency of the signal to be captured (64 MHz) leads to a reduction of unwanted signal noise (Fig. 2.15). Reducing the bandwidth of a factor ∆F increases the The choice of the acquisition plane has also been analysed. As illustrated in Sec. 1.1, the hepatic artery sprouts from the celiac trunk and then reaches the liver, situated on the right side of human body. In a simplified vision we can consider the hepatic artery like a straight vessel going from the centre of the body to the right. Actually, in most cases it goes slightly down and then up again when entering the porta hepatis. Even if axial and sagittal images were tested, we have finally chosen to acquire the images on the coronal plane for a better visualisation of the hepatic artery.

In order to respect clinical times, we reduced the number of slices. The choice of 30 slices implied the exclusion of a large part of the body from the FOV, but this did not lead to important artefacts around the region of interest (the hepatic hilum). In our case, the choice of the phase encoding direction giving less artefacts was the right to left direction.

2D phase contrast optimisation

The principal advantage of 2D imaging is the strongly reduced acquisition time, that allows for breath hold acquisition and much more reliable quantifications. On the other side, it inevitably furnished less data than an entire 3D image: only one segment of the hepatic artery has to be chosen.

Two dimensional acquisitions were tested on each of the three standard directions. The correct placement of the imaging plane is not trivial. The hepatic artery was not visible on any of the reference images acquired during our study without contrast agent injection, including T1 and T2 weighted imaging. Even in clinical routine, a rapid and efficient choice of the optimal plane can be critical. We first realised five 2D images on the coronal plane with a 5 mm slice thickness every 10 mm, but this resulted in a too high loss of time and a partial volume effects compromised imaging. The final choice was to acquire images on the sagittal plane, increasing the chances to easily include the hepatic artery in an imaging plane perpendicular to it.

The choice of sagittal imaging also allowed us to increase the signal to noise ratio, through the choice of a larger flip angle and slice thickness than in 3D. Indeed, if the vessel is perpendicular to the imaging plane, the risk of signal saturation is strongly reduced. The protons that are excited by the θ pulse quickly exit the field of view and new protons with not yet flipped magnetisation vector come in. Moreover, partial volume effects are avoided if the slice thickness increases in the same direction of the vessel. The selection of the imaging plane can be more easily realised in clinical routine if based on the image acquired at early arterial phase, when the hepatic artery enhancement is the highest and its visibility greatly ameliorated. In order to facilitate the finding of the hepatic artery, the sprouting of the celiac trunk from aorta should first be identified and then followed until its bifurcation into splenic artery and common hepatic artery.

Drawbacks of this choice of imaging plane include the complex recognition of the vessel of interest on a slice that is not commonly visualised in clinical routine in the absence of neighbouring slices. This hindrance can be addressed thanks to the comparison and registration, available in most imaging softwares, with any other protocol image reconstructed on the sagittal plane. In particular the comparison to the sagittal reconstruction of the 3D PC sequence appears to be more useful for the comparison of vascular structures, due to their very similar contrasts.

During these acquisitions, the test object was placed under the back of the patient, approximately behind the liver. Saline solution was injected at a flow rate of 0.5 ml/s.

Being the diameter of the tube supposed to be equal to 3 mm, the liquid should flow at a velocity of 7.07 cm/s. The parameters of the 2D and 3D phase contrast optimised sequences we propose are given in Table 2.2.

Results on test object

The employed test object and an example of the acquired images are shown in Fig. 2.16.

We recall that the aim of this part of the study was to validate the post processing of the

V all = V 2 x + V 2 y + V 2 z .
The average of the measurements of absolute velocity in different sections of the tube is 7±0.4 cm/s.

Results on volunteers

The optimised 3D and 2D phase contrast MRI sequence, whose parameters are given in Table 2.2, were acquired on a volunteer and the obtained images are shown in Fig. 2.17 (Fig. 2.17), many more vessels than in the single 2D slice are visible. Blood velocity was measured on aorta, vena cava, and celiac trunk for validation purposes. The velocity in the common hepatic artery and in the left hepatic artery are also given in Table 2.3. Concerning the 2D sequence, in Fig. 2.18 the three images giving a contrast proportional to the velocity on one of the three standard directions are shown first. Below, the magnitude-weighted phase image produced by the GE MR system is compared to the magnitude image and to the reconstructed velocity map. The latter was calculated thanks to the equation (2.1) and its grey value corresponds directly to the velocity in mm/s. The vessels visible on the sagittal 2D image are the hepatic artery, the portal vein, recognisable because it enters in the hepatic hilum next to the hepatic artery, and the renal artery, whose contrast is highly enhanced thanks to the great blood amount of blood flowing to the kidneys. The measurements performed on such images on the mentioned vessels and on the test object are shown in Table 2.4.

Results on patients

Due to the encouraging results obtained on the volunteers, the optimised 2D and 3D A more complicated case is proposed in Fig. 2.20. The celiac trunk and the beginning of the common hepatic artery have a distinguishable contrast, but even the first bifurcation to the gastroduodenal artery is non visible, nor the proper, left and right hepatic arteries.

The measured values on the principal vessels of all the considered patients are given in Table 2.5 for the 3D images and, where available, the 2D ones. indicate that the vessel could not be recognised or that it was not included in the image field of view. For example, in the case of patient 2, illustrated in Fig. 2.20, proper, left and right hepatic artery were not visible.

Validation

Measurements validation was performed thanks to the injection of saline solution in the tube placed under the back of the volunteer. Moreover, it was possible to test ultrasonography measurement of blood velocity on one of the volunteers. In this case, only blood velocities in descending aorta and vena cava were compared, as shown in Table 2.6.

It can be noticed that PC MRI gives good performance for the measurement of blood velocity in aorta, and overestimates velocity in the vena cava.

Discussion and conclusion

Two MRI contrast phase sequences optimised for the measure of blood velocity in human hepatic artery are proposed in this chapter. The sequence has been tested on a test object, on healthy volunteers and is now realised on patients at CEM in Rennes. The measurements of velocity on the images of the test object were validated against ground truth. This let us verify the computation of the velocity map. The difference between real (6 cm/s) and measured velocities (7±0.4 cm/s) is indeed negligible. Blood velocity measurements on volunteers were also validated thanks to the comparison to the velocity of the saline solution injected in the test object placed under the back of the volunteer.

Due to the coherent results obtained, the sequence was applied on HCC patients. The results finally obtained on the 8 patients appear to be coherent. As already told in Chapter 1, in the literature blood velocity in common hepatic artery varies approximately between 10 and 50 cm/s. In particular, [START_REF] Dyvorne | Quantification of hepatic blood flow using a highresolution phase-contrast MRI sequence with compressed sensing acceleration[END_REF] assess the velocity of hepatic artery to 11.6 or 13.7 cm/s, based on two different techniques of phase contrast MRI. The decreasing of the average blood velocity from celiac trunk to the common, proper, right and left hepatic artery is also coherent with the decreasing of the relative diameter of such vessels.

The measurements acquired on portal vein also appear to be coherent with the literature. For example, [START_REF] Carlisle | Estimation of total hepatic blood flow by duplex ultrasound[END_REF] assess the average blood velocity in portal vein to 14.2 cm/s, whereas [START_REF] Dyvorne | Quantification of hepatic blood flow using a highresolution phase-contrast MRI sequence with compressed sensing acceleration[END_REF] propose a value equal to 10 cm/s, based on phase contrast MRI. The values measured on PC MRI by [START_REF] Dyvorne | Quantification of hepatic blood flow using a highresolution phase-contrast MRI sequence with compressed sensing acceleration[END_REF] in hepatic artery, and in particular in portal vein are reasonably close to the values we propose in Table 2.5.

Studies addressing velocity quantification in aorta and portal vein show an average velocity in aorta to 10 cm/s [Gabe et al., 1969, Segadal and[START_REF] Segadal | Blood velocity distribution in the human ascending aorta[END_REF] and in vena cava to 14 cm/s [START_REF] Gabe | Measurement of Instantaneous Blood Flow Velocity and Pressure in Conscious Man with a Catheter-Tip Velocity Probe[END_REF]. In such vessels, the values we have measured with MRI are less consistent with the literature and with the relative values measured with Doppler US. It should be said though, that the sequences we propose have been optimised for small vessels like hepatic artery. Consequently, the velocity encoding value may be unadapted to estimate of blood flow in larger vessels like vena cava and aorta, where it is known to follow a different velocity profile. Moreover, the field of view is centred on the hepatic artery and artefacts increase while leaving the centre of the FOV. Imaging plane and phase encoding directions were also chosen in order to avoid artefacts on the hepatic artery, which may thus alter the imaging contrast elsewhere.

In Table 2.5 it can be noticed that velocity values in the 2D image are generally higher than those in the 3D ones. The average difference is higher concerning values in the hepatic artery (11 cm/s) and decreases in portal vein (4 cm/s) and renal artery (3 cm/s). This is not surprising, the 3D sequence cannot be acquired during a breath hold, and image quality is strongly reduced due to breathing movements. The smaller error in the renal arteries can be related to the higher distance from diaphragm of kidneys than liver, and the consequent decreasing of movements related to breath. Hepatic artery and portal vein are very close and the velocity differences cannot be here related to the distance from diaphragm. The smaller error in velocity measurements on portal vein can be related to its bigger diameter, that may reduce partial volume effects related to breathing.

An additional validation process will nevertheless be needed; a regular comparison between phase contrast and Doppler US is planned for future patients that will undergo such PC MRI. Even if this technique still needs to be optimised, a method for the validation of numerical simulations of blood flow in patient-specific hepatic artery was set up and blood velocity was measured on 8 HCC patients so far.

Chapter 3

At-risk patients identification

Among the more than 200 diagnostic angiographies realised at the Eugène Marquis regional Centre for the fight against cancer that we could analyse, in only 97 cases the patient could actually benefit from the treatment, due to multiple and various reasons.

In this chapter medical image processing is applied to an earlier identification of possible elements that might prevent the patient to benefit from the Theraspheres ® treatment. An earlier exclusion from the protocol would avoid the performance of an invasive and costly surgical procedure. Furthermore, the earlier an alternative treatment is planned, the better will be the life expectancy of the patient. In order to deliver to the tumour a radiations dose high enough, healthy parenchyma is often irradiated as well, usually under the limit threshold of 120 Gray. Effects on liver parenchyma include fatigue and more rarely, abdominal pain and slight and temporary alterations of liver function. Serious liver decompensation (dysfunction of hepatic activity)

is very rare and risk factors are related to pre-treatment gallbladder disease, portal vein thrombosis and a dose to the healthy parenchyma higher than 120 Gray.

Extra hepatic deposition does not only occur, as already mentioned, to the lungs, but microspheres might flow through the gastroduodenal or gastric arteries (cf. Nevertheless, hepatopulmonary shunt is a major contraindication due to the impossibility to treat it directly: the tumour induced capillaries that open a direct way from arteries to veins cannot be precisely identified, and therefore no surgery can lead to their occlusion. The following section is devoted to an extensive description of hepatopulmonary shunt and its possible remedies.

Hepatopulmonary shunt

When a tumour grows into the hepatic tissue, it soon activates an angiogenic process in order to have a direct access to nutrients. Such malignant new vessels can create a short circuit between different vascular networks that are normally supposed to be connected only through the endothelial tissue. Arterioportal and arteriovenous shunts can occur in liver. The latter can also be induced by a portal vein thrombosis or by the arterial embolisation realised during a possible previous chemoembolisation treatment. Such vascular modifications induce pressure variations among the vascular network, and can lead to the development of shunts [START_REF] Ngan | Arteriovenous Shunting in Hepatocellular Carcinoma : Its Prevalence and Clinical Significance[END_REF]. When the radioactive microspheres are injected through the hepatic artery, they normally get stuck into the small capillaries and irradiate the surrounding tissues. In the presence of an arteriovenous hepatic shunt (Fig. 3.1), some of the microspheres may flow directly into the hepatic veins, reach the inferior vena cava and then the lungs. This organ can tolerate no more than 30 Gy for show that these thresholds will be passed, the patient must be excluded from the protocol.

Currently, when an arteriovenous shunt is identified, physician can actually undertake alternative procedures that might lead to a reduction of the shunt. First, the amount of the dose to be injected can be reconsidered. A compromise has to be reached between the radiations limit to the lungs and a still effective dose to the tumour. Reducing the dose amount can ensure a tolerable radiations concentration to the lungs, but on the other hand the dose to the tumour can become insignificant to justify the treatment. Another option that can be considered if the origin of the shunt is identified, is to treat the shunt itself. An anticancer drug like Sorafenib includes anti-angiogenic agents that may indirectly reduce the lung shunt fraction. In other cases, the distribution of the injected therapy can be controlled by either choosing a more convenient injection point (and thus performing a second diagnostic angiography), or occluding a particular artery.

State of the art about shunt quantification

Currently, every hospital centre performing HCC treatment with Theraspheres ® or SIR-Spheres ® relies upon the planar scintigraphy, performed after the diagnostic angiography, for the detection and quantification of extra hepatic depositions. In order to compute the shunt fraction, the nuclear physician identifies a different region of interest (ROI) for each lung, for the liver and for the mediastinum, the latter being used to determine the background noise. An example of the selection of different ROI on a 2D lungs scintigraphy is given in Fig. 3.2. The dose activity is measured on each ROI and the percentage present in the lungs is computed. Nevertheless, the possibility to avoid the realisation of the diagnostic angiography is appealing and many research teams have performed interesting studies on this subject. [START_REF] Gaba | Characteristics of Primary and Secondary Hepatic Malignancies Associated with Hepatopulmonary Shunting[END_REF] realised a retrospective study on 141 patients who undertook the first diagnostic angiography, the target being to identify a statistical correlation between some tumour characteristics (tumour type and size, portal vein invasion and arterioportal shunt) and the lung shunt fraction. Different categories of shunt (<10%, 10-20%, > 20%)

were associated to differences in tumour characteristics, but no parameter was shown to establish a linear correlation with the lung shunt fraction.

The existence of a relationship between the opacification of the hepatic veins at early arterial phase after contrast agent injection, and the presence of an arteriovenous shunt has been mentioned in many studies among which [START_REF] Schneider | MRI of the Liver. Imaging Techniques, Contrast Enhancement, Differential Diagnosis[END_REF], Stabile et al.

[2004], [START_REF] Wu | Liver involvement in hereditary hemorrhagic telangiectasia: CT and clinical findings do not correlate in symptomatic patients[END_REF], [START_REF] Mei | Transcatheter Arterial Embolization of Hepatic Arteriovenous Shunts in Patients with Hepatocellular Carcinoma[END_REF].

A link between arteriovenous hepatic shunt and high pulmonary shunt has been shown

by [START_REF] Olorunsola | Imaging Predictors of Elevated Lung Shunt Fraction in Patients Being Considered for Yttrium-90 Radioembolization[END_REF], who performed a retrospective study based on the analysis of multiphase CT-scan images of 117 patients, including 15 diagnosed with pulmonary shunt.

The target was to show a direct correlation between lung shunt fraction and either a tumour thrombus in a hepatic vein, its occlusion due to external tumour compression, or the early opacification of hepatic veins. The first two characteristics were related to a shunt fraction higher than 20 % on no more than 5 patients. The early enhancement of hepatic veins was supposed to be an evidence of tumour induced arteriovenous shunt, and was visually evaluated on the late arterial phase of the multiphase CT-scan. A higher correlation was revealed between the arteriovenous shunt, identified by the opacification of hepatic veins, and a high shunt fraction (>20 %). Nevertheless, this study was based on the analysis of CT images, which is more invasive than MRI and early enhancement of hepatic veins was evaluated only visually by two radiologists, a third one being consulted in case of conflict.

Moreover, the conclusions proposed by [START_REF] Olorunsola | Imaging Predictors of Elevated Lung Shunt Fraction in Patients Being Considered for Yttrium-90 Radioembolization[END_REF] need to be corroborated on a larger cohort of patients diagnosed with hepatopulmonary shunt. Globally, none of these studies appears to propose a valid, non invasive shunt quantification method, able to provide reliable data without the realisation of the scintigraphy post 99m Tc-MAA injection.

Proposed method for the quantification of lung shunt fraction

Hepatopulmonary shunt can occur if cancerous small vessels create a direct link between the hepatic arterial and venous networks. If this happens, blood flow coming from hepatic artery reaches the hepatic veins earlier than through the normal path. This phenomenon is at the basis of the method that we propose in order to identify hepatopulmonary shunt in a non-invasive way. The imaging protocol used at CEM has been retrospectively analysed and the images of patients presenting a hepatopulmonary shunt have been treated. The analysis performed on the available images is described below, and an improvement of the current multi-arterial MRI sequence is proposed in section 3.7.2.

In order to verify our hypothesis, a monitoring of contrast product distribution was performed through a semi-automatic drawing of ROIs on the hepatic veins at different phases of contrast agent distribution. In Fig. 3.3 is presented a case of a remarkable enhancement of hepatic veins during the arterial phase. It can be noticed that during this phase, the hypervascular tumour is enhanced, as well as aorta, and that vena cava is not reached by the contrast agent yet. We compared the lung shunt fraction to the measure of the different contrast uptake in the same region of interest, before and right after the injection.

In order to improve the quality of the image analysis, we compared the contrast uptake of the hepatic veins not only over time, but also at different positions. Indeed, a reflux of contrast agent towards the inferior vena cava, and even in the hepatic veins, is not rare in patients affected by tumour diseases. The contrast media injected into the arm vein reaches the heart approximately 10 seconds after the injection. If the patient suffers from cardiac diseases, the contrast agent may flow down to the inferior vena cava instead of entering the right atrium. The presence of contrast agent into the vena cava or the hepatic veins could be then misinterpreted as an arteriovenous shunt. This is why we decided to compare different measures taken on two regions of interest placed on the inferior vena cava, the first one before the hepatic veins, and the second one downstream of them. The comparison between these two values can help in the differentiation of an actual pulmonary shunt from a vein reflux disease. If the simple measure of contrast uptake in the hepatic veins does not lead itself to a clear finding concerning the shunt, the second measure can solve this hindrance. Indeed, a contrast agent reflux down to the vena cava would generate a general enhancement of the vessel, independent from the positioning of the hepatic veins:

if the contrast level in the two regions of the vena cava is similar, the reason is probably a venous reflux, otherwise we might be in front of a case of arteriovenous shunt. In Fig. 3.4 a schematic representation with yellow ellipses of the placement of the ROIs in the vena cava and on the hepatic vein is given.

Images analysis

Thanks to the strong partnership with the Eugène Marquis centre, we collected a database made by more than 200 inclusions in the SIRT protocol, meaning that a same patient could be included multiple times in the protocol, for example after an attempt to reduce pulmonary shunt percentage. Among the cases considered, 31 patients had an hepatopulmonary shunt higher than 5% and up to 84%, with an average of 24 ± 20.6%. A control group of 40 patients without pulmonary shunt was also included in the study. Every patient had been submitted to a multiphase CT-scan examination, presenting at least three kinds of images: without intravenous injection, at the arterial phase, and at the portal phase. An MR examination is not currently always included in the protocol, even if its use is constantly increasing thanks to continuous research advances. The possibility to perform multi-arterial acquisitions with MRI is fundamental in particular for the target of this study. Therefore, all the patients to whom multi-arterial MRI acquisitions were previously performed, within the Theraspheres ® framework, were retrospectively analysed as well. The examination consists of one image acquired before the injection, three during the arterial phase, one at the portal phase, another one three minutes after the injection and the last at the late phase.

Analysis of the multiphase CT-scan acquisitions

Pixel by pixel follow up of the contrast uptake was not possible due to important differences in acquisitions parameters that complicated the registration of the same patient's images. Changes in spatial resolution, field of view and number of slices are motivated by the aim to reduce the radiations dose delivered to the patient. Among the 71 cases included in the protocol, 12 were removed (among which 2 with a pulmonary shunt), due to an inconsistency in the acquisitions parameters of images imported from other hospitals, invalidating the coherence of measurements.

We decided to draw separate regions of interest on every image. The ROIs were placed not only on the vena cava and hepatic veins, but also on healthy parenchyma, aorta and vertebral muscle (cf. Fig. 3.3) in order to perform a statistical analysis on a larger number of data and to always have reference values for normalisation.

The placement of the ROI on the hepatic vein and on the vena cava presented significant difficulties due to their poor visibility during the arterial phase and before the injection.

Moreover, advanced stage HCCs happen to compress or even invade hepatic veins or vena cava, making it more challenging to distinguish them and in some cases not enough clear to allow any measurement (cf. Fig. 3.5 for the vena cava and Fig. 3.13 for the hepatic veins). Concerning in particular the two ROIs on the vena cava, it is important to draw The contrast uptake ratio between the non injected phase and the portal phase ∆C n p , as well as between the arterial and portal phases ∆C a p were evaluated as well in order to quantify the contrast kinetic behaviour:

∆C a p = vp -va va and ∆C n p = vp -vn vn , (3.2)
where vp is the normalised value measured at the portal phase. The ratio between these values was also computed as

∆C p = ∆C a p ∆C n p . (3.3)
The idea justifying the analysis of this value is the need to describe the contrast uptake not only until the arterial phase, but to provide a global vision of the contrast distribution during the three phases. A value of ∆C p close to zero is obtained if the contrast uptake from the arterial phase to the portal phase (∆C a p ) is small and the one from the non injected phase to the portal phase (∆C n p ) is high. This would mean that most of the contrast agent reaches hepatic veins during the arterial phase, suggesting the presence of an important shunt. Otherwise, if the enhancement between the arterial phase and the portal one (∆C a p ) is much bigger than ∆C n p , most of the contrast agent reaches the hepatic veins during the portal phase, which is normal, and the shunt percentage may be negligible.

Finally, the last measures concerned the difference between the average grey value of the two ROIs placed on the vena cava (cf. Fig. 3.6, on the bottom right and top right images), the first one placed before the arrival of the hepatic veins (v bHV ) and the other after it (v aHV ):

∆C V C = vaHV -vbHV vbHV .
(3.4)

Analysis of the multi-arterial MR images

A proper statistical study could not be performed on the MRI acquisitions due to the low number of patients presenting a hepatopulmonary shunt higher than 5 %. Among the 18 patients who underwent an MR examination, two had important pulmonary shunt fractions equal to 15.7% in one case and to 22.7% in the other case.

MR images were nevertheless analysed due to the strong importance we think this kind of imaging will have in the future. MRI provides a better image contrast than CTscan, and therefore an easier visualisation of the vascular networks, despite a lower spatial resolution. The acquisition of three image sets during the arterial phase is crucial for a correct monitoring of the kinetics of the contrast agent in the vessels, and thus of the blood.

Such finer temporal resolution allows as well to easily differentiate a venous reflux from an arteriovenous shunt, thanks to the simple analysis of the contrast dynamic in the veins.

A single acquisition during the arterial phase is a strongly limiting factor of CT-scan, and our study retrospectively showed that the contrast agent distribution was captured at slightly different moments for different patients, leading to possible inconsistencies in measurements.

An example of the series of MRI acquisitions during the circulation of contrast agent is given in Fig. 3.7. The MRI sequence used for the follow up of contrast agent distribution is the LAVA (Liver Acquisition with Volume Acceleration) sequence, already described in section 2.2.2: it is a T1-weighted gradient echo sequence currently performed at the injection time (Fig. 3.7b), 9 seconds later (Fig. 3.7c) and 18 seconds later (Fig. 3.7d).

Only these images were analysed because the acquisition before the injection (Fig. 3.7a)

and the ones at portal (Fig. 3.7e) and late phases are realised using the slightly different sequence LAVA-Flex, which has a longer acquisition time but gives better contrast and spatial resolution. A quantitative comparison between images acquired with different sequences being not reliable, the first image of the arterial phase was considered as a reference, due to the presence of contrast agent only in the aorta (Fig. 3.7b). Thanks to the multiple acquisitions during the arterial phase, a curve representing the hepatic veins enhancement in time can be drawn.

The ROI placement was realised following the same criteria as in the analysis of the CT-scan images. Normalisation was performed measuring the grey value of a ROI placed on an empty region (in air) in order to measure the noise of the image.

An example of the drawing of ROIs on the three hepatic veins (right, middle and left)

is given in Fig. 3.8. The ROIs were semi-automatically drawn thanks to the definition of macros in the Fiji software and their placement was validated by an experienced radiologist.

Statistical analysis of the image based data

For a patient population P, a same statistical analysis was performed on each of the 3 parameters previously defined (∆C n a , ∆C p and ∆C V C ). This analysis involves:

• a statistical Student's t-test, to detect if there is a significant difference or not be-tween the parameter average values, computed on the subpopulation P 1 (including the patients with a shunt value superior or equal to a given threshold) and on the complementary population P 2 ( P = P 1 ∪ P 2 ),

• a linear regression analysis between the parameter (the predictive variable) and the shunt percentage (the variable to be predicted),

• the identification of a supervised classifier, trained to automatically decide, based on the computed parameter value, if patient's shunt value is superior or not to the chosen threshold.

K-nearest neighbours classification algorithm (with k =1) was used. A cross-validation procedure was used to estimate the classifier performance, by randomly partitioning the population P into 10 subsets P 1 , . . . , P 10 , each of them including approximately the same number of patients. The classifier was trained 10 times on the subsets

P i = 10 j=1 j =i P j , i = 1, . . . , 10.
The classifier learned on P i was applied to the test population, and the performances (specificity, sensitivity, classification error) were computed on all the results obtained on the 10 tested subpopulations.

The partitioning of the population was also used for the correlation analysis: the correlation coefficient between the parameter and the shunt percentage was estimated 10 times, respectively on P i , i = 1, . . . , 10, in order to obtain an evaluation of the corresponding estimation dispersion. An average of such 10 estimations is then computed.

Concerning the analysis of CT-scan, 59 patients were included in the study, 29 presenting a shunt percentage higher than 5%, and 30 without a significant shunt. Statistical analysis was performed using the softwares Weka (University of Waikato, New Zealand)

and Matlab 2015b (The MathWorks Inc., Massachusetts, USA.)

It should be noticed that the measured parameters correspon to an average, computed on a number N p of pixels in a region of interest. The standard deviation for a given parameter can be approximated by σ √

Np

, where σ is the standard deviation computed on the N p values.

Student's t-test

In the case of SIR-Spheres ® , a strict contraindication to the treatment is a shunt percentage higher than 20%. The limit threshold for the Theraspheres ® treatment is rather related to the radiation dose delivered to the tumour. Nevertheless, a 20% shunt often contraindicates both treatments. Therefore, the shunt percentage threshold for patients classification was defined as 20%.

Since the patients included in our study had undertook a diagnostic angiography in view of a possible Theraspheres ® treatment, the ground truth value of shunt percentage was available from the analysis of the post angiography planar scintigraphy. Among them, 46 patients had a shunt percentage lower than 20%, and 13 higher than 20%.

The procedure described below was applied for each of the three measured parameters ∆C n a , ∆C p and ∆C V C . The gathered data were divided into the two classes C 1 (< 20%) and C 2 (>20%). The Student's t-test (significance level equal to 0.05) was applied in order to confirm or not the sensitivity of the parameter statistical average, to the class factor (C 1 or C 2 ). Every class was preliminary verified to follow a Gaussian distribution using a one-sample Kolmogorov-Smirnov test. Different Student's t-tests should be used if the variances of the two populations to be studied are equal or not. The equality of populations variances was thus investigated using a F-test and the appropriate Student's t-test was applied to each variable.

Classifier validation and correlation measure uncertainty

Cross-validation and the k-nearest neighbours method (k=1) are used in order to measure the classifiers performance. The population has been randomly partitioned into 10 subsamples P 1 , . . . , P 10 . The following procedure was applied ten times: one subsample is removed and the remaining set is used to decide to which category every elements of the removed subsample should belong. The chosen category corresponds to the one that prevails among the k closest neighbours of the considered element. The value of k is empirically chosen as the one providing the best results, and is defined to 1 in our study. For every classified element, the computed class is compared to ground truth. At the end of the last iteration, the sensitivity, specificity and classification error are computed.

Concerning the computation of a correlation coefficient, it was computed using the same partitioning procedure presented above. It is computed 10 times on 9 subsamples (at 3.6. RESULTS every iteration, one subsample is removed from the total population). The 10 results are averaged resulting in a final correlation estimation. The dispersion is evaluated, using the bootstrap approach.

Results

Computed tomography imaging

The average and standard deviation values of the three measures contrast values for each of the two populations is shown in Table 3.1 Table 3.1: For each of the two populations, this table presents mean and standard deviation of ∆C n a (the hepatic veins enhancement between non-injected phase and the arterial phase), ∆C p (the ratio between the non injected-portal phase enhancement and arterial-portal phases enhancement) and ∆C V C (the difference between the contrast uptake on the vena cava before and after the hepatic veins at the arterial phase)

∆C n a ∆C p ∆C V C
Patients with shunt < 20% 0.5 ± 0.4 0.5 ± 0.2 0.05 ± 0.3

Patients with shunt > 20% 1.8 ± 1 0.05 ± 0.2 0.4 ± 0.7

The first step before the application of statistical tool is the visualisation of the gathered data.

Thanks to the measurement uncertainty study realised on the average grey value measured in each ROI, every measurement will be visualised through its average and standard error. The three principal measured values were plotted, as a function of the shunt percentage. The most meaningful graph is given in Fig. 3.9, which represents the values of the hepatic veins enhancement at the arterial time (∆C n a ) versus the shunt fraction and appears to confirm the hypothesis of an existing correlation. The best linear fit computed with a linear regression is also shown. The same graphs are also shown for the other two considered values in Fig. 3.10 and Fig. 3.11.

The linear regression results concerning all the values measured on the CT-scan data are shown in Table 3.2: the hepatic veins enhancement between non-injected phase and the arterial phase ∆C n a , the ratio between the non injected-portal phase enhancement and arterial-portal phases enhancement ∆C p and the difference between the contrast uptake on the vena cava before and after the hepatic veins at the arterial phase ∆C V C . The table gives the correlation coefficients as well as the equations of the linear regression lines which 3.6. RESULTS best fit the data of the shunt and of the measured enhancements. ) and the difference between the contrast in the vena cava before and after the hepatic veins (∆C V C ). On the contrary, population variances were estimated as equal in the remaining case, and the consequent Student's t-test was applied to each different value. These results are summarised in Table 3.3, together with the p-value obtained by the realisation of the t-Student's tests, a p-value smaller than 0.05 generally revealing good statistical significance. Predicted class: Shunt < 20% 5 41

Magnetic resonance imaging

The better visibility of vascular network in MR images allowed us to easily identify all the three hepatic veins (left, median and right) on each of the three images acquired during the arterial phase, on 9 patients. This important information has been used in order to realise enhancement curves, describing the contrast uptake of each of the hepatic veins over time. The plot of a different curve for each hepatic vein allows identifying the possible presence of a shunt in case only one of them presents an early enhancement. In particular, the enhancement curves concerning the first patient presented in Fig. 3.12 display a clear enhancement of the median hepatic vein, not visible on the other veins.

Discussion and perspectives

Globally, the analysis of the CT images of 59 patient whose shunt percentage was known, let us confirm, by the regression shown in Fig. 3.9 and the Student's t-test, the hypothesis concerning a relationship between the pulmonary shunt percentage and the early enhancement of hepatic veins.

Nevertheless, results do not allow concluding the existence of a clear linear dependency between these two variables. The presence of values unrelated to the best approximation of the linear regression line can not be neglected. The reliability of the hepatic veins enhancement as a parameter determining the shunt fraction is still limited. For the moment, a quantification of the early opacification of hepatic veins does not appear to be enough reliable to estimate the lung shunt fraction and thus to replace the diagnostic angiography and scintigraphy. Concerning the analysis of MR images, a statistical correlation has not been investigated due to the low number of patients with a shunt. Nevertheless, a better visualisation of the vascular network, as well as a better temporal resolution could lead to a more adequate estimation of the shunt percentage.

Limitations of the current study

Some important restriction to this study on CT and MR images must be considered, in addition to the low number of patients with a pulmonary shunt percentage higher than 20%, compared to those with a lower shunt percentage.

First, the presence of a venous reflux can induce an early enhancement of hepatic veins not related to a shunt, but possibly to the large blood amount flowing to the kidneys, that can help the diffusion of contrast agent early in the vena cava. The distinction between this phenomenon and a shunt is critical, but impossible to perform if only one image is acquired during the arterial phase. The importance of MRI is once again crucial. Indeed, only a tracking of contrast agent over time can help identifying the origin of an early enhancement of hepatic veins. The possible distinction of the three hepatic veins is also fundamental: a venous reflux concerns usually all the three hepatic veins, whereas a pulmonary shunt can be identified by the different uptake in the three hepatic veins.

Second, the presence of massive or infiltrating tumours can highly complicate the finding of hepatic veins and affect the measurements. An example of such kind of obstacle is given in Fig. 3.13, which shows a CT-scan acquired during the arterial phase, where the presence of an extensive tumour prevent the localisation of the mean hepatic vein. A solution might rely in the decision to realise, for the patient concerned, a more detailed image analysis aiming to perform an advanced description of the vessels geometry, for example through dynamic MRI. Third, MRI analysis has shown to be more reliable that CT, even if some incoherences exist, especially concerning the different enhancement curves in the absence of a pulmonary shunt (Fig. 3.12c). In order to understand if a similar pattern, showing one of the veins more enhanced than the others, is due to a low image resolution or to an actual shunt, it can be suitable to realise a second diagnostic examination with better sensitivity and resolution, particularly addressed to the shunt detection. In this case, it is important to ensure that the most enhanced vein actually derives from the tumour tissues.

A last possible error concerns patients affected by a high pulmonary shunt, but whose hepatic veins enhancement is low, at least on the only arterial image acquired with CTscan. Such cases demonstrate that a single CT-scan acquisition is not sufficient for a correct detection of shunt. The synchronisation between the injection of contrast agent and the acquisition cannot always ensure to capture its distribution exactly when the contrast agents reaches the hepatic veins.

Concerning the analysis performed on the enhancement of the vena cava during the arterial phase, no significant results have been obtained, and this measure will not be considered in the future, at least for shunt detection.

The statistical significance of this study will have to be investigated through an analysis on a larger patients cohort. Moreover, we envisage to perform an actual statistical study on images acquired with dynamic MRI, which appears to be the most encouraging technique.

In particular, we envisage a revision of the currently employed MRI sequence for the dynamic follow up of contrast distribution, in order to tailor it on shunt detection without losing any of its benefits concerning HCC characterisation. Preliminary investigations have been undertook and are described in the following section.

Proposed MRI sequence for the detection of hepatopulmonary shunt

The flexible and uninvasive characters of MRI are important advantages even in the present context. The restriction to axial imaging in the case of CT is likely to have an influence on the preponderance of this kind of acquisition even in other imaging techniques. MR images are often acquired on the axial plane for simple reasons related to clinical practice, influenced by the familiarity of physicians to the axial vision of human body. Nevertheless, flexibility in the choice of the acquisition plane allows to acquire coronal images, which we consider to be the best choice for the imaging of hepatic veins. These vessels can indeed be approximately considered as three lines coming from the left, central and right sides of the liver and reaching the vena cava just above the liver (cf. Fig. 1.2). Their verticality makes them more visible on coronal images, where a longer segment can be seen on a single slice, allowing for a more comprehensive vision of their geometry.

On the other side, the non-invasive character of MRI allows the operator to realise as many images as necessary, as long as it is compatible with acquisition time. Risks for the patient concern indeed only some rare cases of allergic reactions to contrast agent. If a rapid and frequent monitoring of the enhancement of hepatic veins is achieved, important qualitative and quantitative information about hepatopulmonary shunt can be captured.

This is why we propose the adoption, in clinical routine, of an MRI sequence characterised by multiple arterial acquisitions, as frequent as possible, in the coronal plane, replacing the axial LAVA sequence currently employed. We performed two acquisition tests on a healthy volunteer, without the injection of contrast agent. The LAVA sequence currently realised on the axial plane in the HCC protocol was used as a reference and slightly modified. First, parameters were tuned in order to obtain a very good quality image: in Fig. 3.14a the geometry of hepatic veins is clearly visible, even in the absence of contrast agent. Nevertheless, the resulting acquisition time, equal to 15.6 seconds, can not be accepted in clinical routine. Next, the sequence parameters given in Table 3.5 were tested in order to reduce the acquisition time and increase the possible number of images to be acquired during the arterial phase. In particular, the repetition time was reduced The sequence that we propose to adopt in clinical routine, replacing the currently used axial LAVA sequence, provides images with sufficient SNR and spatial resolution, in a globally reduced acquisition time. We also verified that a coronal acquisition leads to a better visibility of the hepatic veins. Four clustered images could be acquired right after the injection of contrast agent, instead of the three currently realised. This could provide a very accurate description of the kinetics of the contrast agent, increasing the possibilities to detect and quantify hepatopulmonary shunt otherwise than with the diagnostic angiography. Nevertheless, the undertaken investigations are still preliminary and our hypothesis needs to be certified by further analyses. Patient-specific hepatic vascular modelling This chapter opens the section of this thesis dedicated to the mathematical modelling, that will be unfolded in Chapters 4, 5 and 6 and that is outlined in Fig. 1.9. In this chapter, the Box number 7 of the Fig. 1.9 is considered: the simulation of the liver itself, and in particular of its complex vasculature, thanks to a simplified description of blood flow in vessels. A substantially more detailed simulation of blood flow, as well as the simulation of microspheres trajectories in a patient-specific domain, can then be addressed in the following chapters.

The global workflow presented in this chapter, starting from the processing of medical images until the simulation of the arterial vascular tree, is outlined in Fig. 4.1. We propose a procedure where the segmentation of the hepatic arteries is achieved on the 3D conebeam CT angiography, whereas the profiles of liver and tumours are segmented from the CT acquired during the portal phase (Sec. 4.1). The raw data of the segmented images are then treated in order to be imported into the computational model. In particular, the binary image describing the shape of the segmented vessels is submitted to a skeletonisation and converted into a tree structure (Sec. 4.2.1.1). It is finally registered to the portal CT-scan from which the liver is segmented (Sec. 4.2.1.2).

The computational model previously developed in the research team [Kretowski et al., 2003b, Jurczuk et al., 2014] and outlined in Section 4.2.2 is finally used in order to simulate the growth of the arterial tree, which is constrained by the specific patient's data as liver and tumours contours. Thanks to the simulated entire arterial tree, a first simplified

Extraction of patient's data from images

As already mentioned in Chapter 2, the clinical imaging protocol (cf. Fig. 2.8) applied for patients suffering from hepatocellular carcinoma offers different kinds of acquisitions, each one describing particular characteristics concerning the vasculature or the parenchyma.

Unfortunately, there is no available image where satisfactory segmentation of liver shape, tumours and vessels as well can be realised in a common reference system. The 3D conebeam CT angiography is acquired during the arterial phase, when hepatic arteries are characterised by a signal hyperintensity, and vessels down to a diameter of ∼0.5 mm are visible (Fig. 2.6a), but liver contour is hardly detectable. Vice versa, during the portal phase, parenchyma enhancement is maximal but the arteries have been "washed out".

The methods used for arteries segmentation from angiography and liver and tumours seg-mentation from portal CT are presented below.

Arteries segmentation

A wide range of numerical methods exists for vessels segmentation [START_REF] Kirbas | A review of vessel extraction techniques and algorithms[END_REF].

Many of them are generic segmentation methods that have been modified in order to be applied to the vessels case. Nevertheless, these kinds of methods are often inefficient for the segmentation of hepatic vessels, which is noticeably complex [START_REF] Esneault | Liver vessels segmentation using a hybrid geometrical moments/graph cuts method[END_REF].

Hepatic arteries show tortuous geometries, a large number of bifurcations and a wide range of diameters, going from the common hepatic artery (∼5 mm) to the smallest segmentable vessels (∼0.5 mm).

We decided to perform vessels enhancement through an Hessian-based multiscale filter like the 3D Frangi's filter. Methods based on the analysis of the Hessian matrix are indeed shown to provide good segmentations results in most vascular cases, and the multiscale character of Frangi's filter is particularly suitable to the size variability of hepatic vasculature [START_REF] Jimenez-Carretero | 3D Frangi-based lung vessel enhancement filter penalizing airways[END_REF].

The algorithm is based on the idea that, if a given set of voxels contains a tubular structure, then the local second derivative estimated from the voxel grey levels is lower in a direction (i.e. along the direction of the vessel) and higher in the two orthogonal directions (i.e. towards the walls of the vessel). This property can be exploited through an analysis of the estimated eigenvalues (λ 1 , λ 2 , and λ 3 ) of the Hessian matrix. In every voxel such matrix is obtained from second order finite differences, computed on the grey level image previously smoothed with a 3D Gaussian kernel. The size of this kernel is a scale parameter: a larger scale parameter corresponds to a more efficient smoothing. The voxels of an ideal vessel verify:

|λ 1 | ≈ 0 (4.1) |λ 1 | |λ 2 | (4.2) λ 2 ≈ λ 3 . (4.3)
In particular, three geometrical characteristics are calculated from the eigenvalues in order to determine the tubular-like character of a given structure.

R A = |λ 2 | |λ 3 | indicates if it is closer to a plate-like or a line structure. R B = |λ 1 | √ |λ 2 λ 3 |
reflects its divergence from a bloblike form, and S = 3 j=1 λ 2 j is used to distinguish vessels information from background noisy structures (the background noise becomes, after smoothing, a low and approximately constant grey level image, leading to small second order differences values). Frangi's filter computes the product of 3 monotonic nonlinear functions (each of them having a range equal to the real interval [0, 1]) of, respectively R A , R B , and S. The sensitivities of the filter to R A , R B , and S are tuned separately by three parameters denoted α, β, and c. A Frangi's filter output near to zero corresponds to a background region and an output close to 1 reveals a vessel structure. The efficient size of the smoothing parameter, i.e. of the scale parameter, depends on the vessel diameter. Therefore two other parameters have to be defined in this method, in order to take into account the size variability of the vessels: a scale range (depending on the diameters of the smallest and the largest vessel visible in the acquisition), and the scale ratio giving the step used to go through this range.

Vessels segmentation can finally be accomplished by a thresholding of the Frangi's enhanced image, followed by a partitioning of the resulting image into connected regions with a union-find algorithm (function "bwlabeln" in Matlab, where voxels are considered to belong to the same connected components if they are linked by at least one of their 26 neighbours). In order to reduce image noise, the classification into connected components requires a very low threshold, whose value is currently manually chosen based on the enhanced image histogram. The output of the this step corresponds to a 3D image where every connected component is labeled by a different grey level. The connected component corresponding to the arterial tree is then visually identified and extracted. This allows the suppression of every enhanced structure that does not belong to the hepatic artery.

Liver and tumours segmentation

Segmentation of liver shape is a challenging task due mainly to its unclear boundaries with other organs and its irregular shape, in particular in the presence of advanced cirrhosis. Various methods for liver segmentation have been developed during the years [START_REF] Heimann | Comparison and evaluation of methods for liver segmentation from CT datasets[END_REF]. The most important challenge in clinical image processing concerns the computational time needed for providing results in rapid clinical routine. The MESA software, based on deformable models, was used for liver segmentation. The software is presented in [START_REF] Reska | MESA: Complete approach for design and evaluation of segmentation methods using real and simulated tomographic images[END_REF] and was developed at Bialystok University of Technology by research partners. MESA is a platform for design and evaluation of segmentation methods and provides tools for management and visualisation of medical images.

The shape of the liver is extracted with a built-in semi-automatic active surface method [START_REF] Reska | Fast 3D Segmentation of Hepatic Images Combining Region and Boundary Criteria[END_REF] that performs a fast fully 3D segmentation of the volume and requires only small manual corrections. The active surface model is represented as a deformable polygon mesh that evolves under the influence of external forces (related to image features) and internal constraints. As shown in Fig. 4.3, the segmentation process starts with an [START_REF] Reska | Fast 3D Segmentation of Hepatic Images Combining Region and Boundary Criteria[END_REF] initialisation of a spherical surface inside the liver. During the evolution, the surface inflates inside the organ to match its shape. The external image energy enables the model to expand only into regions with the image intensity similar to the intensity inside the initial surface, while an internal curvature energy controls the smoothness of the mesh.

Additionally, the model dynamically modifies its topology using an adaptive subdivision procedure that adds additional polygons to the mesh, which is necessary to maintain high quality of the segmentation.

The evolution is an iterative process that can be stopped automatically by the model when further expansion is not possible, or manually by the operator. After the execution, the segmentation results can be shown as an image stack and can be superimposed to the original 2D data (Fig. 4.16). Final results can require some slight manual corrections, especially near the vessels. Corrections are performed with a 2D paintbrush-like tool integrated in MESA software.

The volume of liver was separately measured on a Siemens radiologist workstation, using the vendor supplied application "Volumétrie" in order to validate the segmentation.

Concerning tumours segmentation, being the malignancies borders less neat, they are segmented with an active contour model on 2D slices, which requires major manual corrections.

Modelling of hepatic arterial tree

The computational model dedicated to the simulation of the liver vasculature allows the user to simulate in the same time the growth of hepatic arteries, portal vein and hepatic veins, after an initialisation phase where initial main vessels must be introduced. In the present work, only the hepatic artery will be considered. The contours of the shape of the liver are also needed in order to spatially delimit the vessels growth. The areas where a malignancy is present have different growth parameters, leading to a denser and more irregular vascularisation. The necrotic tissue can also be simulated by a decrease of the vascularisation until total vanishing of the vessels.

The novelty of the work presented in this section is the patient-specific character of the model, achieved through the introduction of initial data extracted from patient's medical records. The use of pure medical data in computational modelling leads to a certain number of challenges that are addressed in the next section.

Initialisation of vascular growth model

First of all, the binary images obtained from liver and tumours segmentation are merged and exported to a raw file made of 256 × 256 × 256 unsigned 8-bit integers going from 0 to 255, used in the vascular simulation model. As it is presented in the next sections, a more important challenge concerns the format of the extracted binary image of hepatic arteries, as well as the registration between different data gathered from different images.

Definition of a binary tree from segmented vessels

The model considers every vascular tree as a binary tree: given a tree root, it must be possible to go through the tree while knowing, for every node, the coordinates of its father node, if it is a bifurcation or a leaf, and, in the first case, if it has one or two children. In particular, the initialisation file must include the coordinates of all the tree nodes, ordered in the "depth-first" manner as illustrated in Fig. 4.4. Starting from the root, the left child Every bifurcation can have 1 or 2 children except the root that has one child of every node is considered until a leaf is found. Then, the last considered bifurcation is searched and its right child is considered next. In addition to the nodes coordinates, the input file must give the node class (if it is a bifurcation or a leaf, also called cell), the number of the cell, based on the depth-first order, and a binary value defining the presence or not of a left or a right child. The possibility to consider bifurcations with only one child allows to approximate a 3D curvilinear vessel branch as a sequence of straight truncated cones. In Fig. 4.4 the input file associated to the example tree is given as well.

In a first time, tree extraction was achieved through a semi-manual procedure applied to 3D angiography with the EndoSize © software (Therenva SAS, Rennes, France). Since this software was originally conceived for the aortic trunk, the segmentation could not be completely automatic, the structure of the hepatic arterial tree being more complex.

Every vessel of the arterial tree was identified by some critical points manually chosen (Fig. 4.5). The 3D coordinates of every selected point were then automatically saved in a format able to reestablish the binary tree structure of the vessels (tree root, children vessels, tree leaves), needed for the initialisation of the simulation in the vascular model.

We also measured the diameters of every vessel segment with EndoSize. Although the tree structure was easily reconstructed, this procedure presented important drawbacks in terms of time and accuracy. vessels registration on liver must be achieved. The registration procedure is described in the following section.

Registration between CT-scan and 3D angiography

As already mentioned, it is not possible to extract liver, tumour and vessels from a same image. Hepatic arteries were segmented from 3D cone-beam angiography and liver and tumours shapes from portal phase CT image. An example of the CT images of a patient is given in Fig. 4.7.

As this image demonstrates, direct registration from angiography to portal phase CT or MRI is not possible due to the absence of any similar and common information on the two acquisitions: hepatic arteries almost do not appear during the portal phase, and liver contour is hardly visible on angiography. This is why we employ arterial phase image to Fig. 4.7b shows that the biggest vessels of hepatic arterial tree are visible on this kind of acquisition (white arrow) and that they can be segmented.

The registration process developed in order to match information issued from different and Fig. 4.7c shows that the first step is straightforward. Indeed, images are acquired during the same examination, the difference lying in the acquisition time during contrast agent propagation and in the choice of some acquisition parameters. We accomplish this step through a rigid transformation identified automatically by a multimodal registration algorithm based on mutual information metric [START_REF] Pluim | Mutual-information-based registration of medical images: a survey[END_REF]. Registration is verified qualitatively through a checkerboard volume.

Registration of angiography on arterial phase image is more challenging. Indeed, angiography is acquired one to two weeks after CT or MRI. Since liver moves to adapt to body changes, the two images are inevitably different. Furthermore, between the two acquisitions, patients can be subjected to extra-hepatic artery coil embolisation, in order to avoid undesirable radiation depositions. Namely, a coil is guided to the artery through a catheter and deployed in the chosen vessel with the aim of occluding it. This inevitably induces some artefacts in the resulting cone-beam CT, which is the case for one of the considered patients (Fig. 4.7a).

The registration is based on the geometry of the main arteries. We draw on both images a region of interest in order to select the biggest hepatic arteries, which are clearly visible even on arterial phase images. Vessels enhancement on both arterial phase image and angiography is then realised by applying a 3D multiscale Frangi's filter [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF][START_REF] Jimenez-Carretero | 3D Frangi-based lung vessel enhancement filter penalizing airways[END_REF]. For the moment, the algorithm has not been tested on MR images. Next to the Frangi's filter, a grayscale threshold is applied to arterial phase CT in order to completely remove every non vascular structure. Last, a 3D erosion is applied on angiography in order to simplify the geometry and remove unnecessary information absent on arterial CT. Images are finally registered through an affine transformation (translation, rotation and scaling). The transformation matrix is automatically calculated by a multimodal registration algorithm, based on mutual information metric [START_REF] Pluim | Mutual-information-based registration of medical images: a survey[END_REF]. The algorithm optimizes images similarity through an iterative evolutionary algorithm [START_REF] Styner | Parametric estimate of intensity inhomogeneities applied to MRI[END_REF] for which we identified a set of suitable parameters. Images resolution can be given to the algorithm in order to identify scaling coefficients, ameliorate results and reduce computational time.

Registration is validated through a visual comparison and superimposition of Frangi's enhanced angiography and arterial phase CT. For further validation, we apply the found transformation to original angiography and visually compare it to original arterial CT-scan.

The identified series of geometrical transformations is finally applied to the list of points describing the arteries, obtained with the procedure described in section 4.2.1.1. Transformation 2 (cf. Fig. 4.8) registers those points to arterial-phase CT, and next transformation 1 is applied to register them to portal-phase CT, where tumours and liver are segmented.

Vascular growth simulation

A detailed description of the vascular growth model, developed in the research team, is

given in [Kretowski et al., 2003b, Jurczuk et al., 2014]. A vascular tree is made of vessels that can divide creating bifurcations (Fig. 4.9a). Each vessel segment (the part of the [START_REF] Zamir | Branching characteristics of human coronary arteries[END_REF], it is assumed that the vascular tree forms a binary tree.

The simulation is initialised with a limited vascular network. First, the size of initial vessels and liver are proportionally rescaled based on a chosen factor included in the interval ]0, 1]. Further development of the vascular trees is modelled as an analogy to a hyperplasia process, namely a growing of tissue amount. The increasing needs of the growing tissue induce the development of the vascular network which is responsible for blood delivery.

After the growth process, vessels and liver reach their original, adult size.

Tissue is represented by a set of macro-cells that are regularly (but randomly) distributed inside the predefined 3D shape (Fig. 4.9b). A macro-cell corresponds to a small part of the tissue, with a given size. The need of each macro-cell for blood is defined by its class, which determines most of its structural and functional properties like its size and birth/death probability, as well as physiological features like blood pressure or blood flow rate. For example, tumour macro-cells have a higher birth probability than the healthy ones, and their spatial density is higher. Moreover, different classes allow various pathological changes to be investigated. Some of the macro-cell features are described by statistical distributions in order to introduce a natural variability. The geometry of smallest vessels like capillaries is neglected. They are treated as to be included in the tissue macro-cells.

As a result, anastomoses (mutual vessel connections), occurring mainly among vessels with very small radii, are not simulated. Thus, the procedure of sprouting new vessels towards the new macro-cell is performed.

Several candidate bifurcations, constrained by a maximal distance from the macro-cell, are created. Based on them, temporary vascular structures appear, as illustrated in Fig. 4.11.

This process can be seen as a competition, since only the vessel that locally minimises the sum of vascular volumes is finally chosen as the permanent one. Afterwards, the characteristics of the vessels, as blood flow and pressure, are recomputed to fulfil all the required physical constraints.

Two physical constraints, represented in Fig. 4.12, are in effect at every bifurcation. The first one deals with the decreasing vessel radii in the vascular tree, creating a relationship between the radius of a father vessel r and the radii of its two descendants r l and r r [Kamiya , 1972]. Second, the matter conservation law must be observed for blood flow at each bifurcation: the physical quantities of blood entering (Q) and leaving a bifurcation

(Q l + Q r )
have to be equal. Moreover, blood flow is modelled in a simplified manner. It is considered as a Newtonian (with constant viscosity µ), incompressible fluid, and its flow is governed by the Poiseuille's law, which relates the pressure difference between the two extremities of a vessel segment (∆P ) to its blood flow rate (Q), length (l) and radius (r) (cf. Fig. 4.12). The method proposed in [Kretowski et al., 2003a] to ensure that these three laws are verified at each bifurcation (the bifurcation's consistency) is detailed in Fig. 4.13.

This method can be extended to ensure the consistency of a structure in which there are two consistent subtrees. In this case, the first step still consists in calculating the pressures blood flow [START_REF] Folkman | Tumor angiogenesis: therapeutic implications[END_REF]. The capillary network is denser even inside the lobules, the elementary components of the liver. In the model, tumour macro-cells and vasculature are defined through three principal characteristics. First, blood flow at the entrance of tumour macro-cells was defined as 3 times higher than in the healthy tissue. This ratio is usually considered to be at least 2 in the literature, but can be much higher [START_REF] Andrews | Modulation of liver tumor blood flow with hepatic arterial epinephrine: a SPECT study[END_REF].

Second, the ratio between birth and death probability of tumour macro-cells was higher than in healthy macro-cells. Last, a denser vascular network was induced through an increase of the spatial density of macro-cells to be irrigated. This is done by authorising a smaller distance between a newly born macro-cell and the already existing ones [Kretowski et al., 2003b]. In this way, the number of macro-cells in a same region is higher, and the size of macro-cells is thus implicitly reduced.

A preliminary simulation of microspheres distribution

We improved the model in order to allow for a better simulation of liver functionalities, in the aim of SIRT modelling. A preliminary simulation of the choice of microspheres best injection point in the arterial tree can be performed. It is temporarily assumed that microspheres flow in every vessel situated downstream from the injection point, proportionally

to the blood flow given by Poiseuille's law. Indeed, an intrinsic property of the vascular model is the notion of blood flow, given by Poiseuile's law, in every vessel and at the input of every macro-cell of the vascular tree. This property easily allows us to compute dose distribution everywhere in the arterial tree. Namely, starting from any injection point in the tree, at every bifurcation met, the dose will be distributed proportionally to the blood flows of the two child vessels. Finally we can calculate, based on any chosen injection point, the percentage of the injected microspheres which reached each healthy and tumour macro-cell. We developed an interactive software that allows the user to navigate through the tree bifurcations, while an interface continuously shows the percentage of irradiated tumour and healthy tissue, as well as the percentage of the microspheres distributed in tumour and healthy tissue.

Results

The complete procedure proposed has been applied to the CT volumes of 2 patients so far.

Vessels segmentation and binary tree extraction have been applied to the 3D angiographies of 2 more patients. Parameters of the anonymised acquisitions used for liver and vessels segmentation from the first 2 patients are listed in Table 4.1. Not surprisingly, the binary tree structure deduced by the analysis of the vessels segmented through the Frangi's filter is more accurate and provides a more detailed initial tree than the one obtained with EndoSize. Concerning the choice of the Frangi's filter

RESULTS

parameters, empirical optimisation was achieved based on the first two patients whose angiographies were analysed. The defined set of parameters appeared to be suitable for every other patient and is given in Table 4.2. For the registration of angiography to arterial phase CT (2 • Registration in Fig. 4.8), a

Frangi's filter vessel enhancement was needed on both images. Frangi's filter parameters for segmentation on 3D angiography have already been given in Table 4.

Vascular growth results for 2 patients

Once all the needed patients information have been collected in a coherent way, the vascular model can be launched based on the patient-specific data. Initially, the binary tree extracted with Endosize have been used. The initial parameters that were used for the simulation of the vascular network are given in Table 4.4. Since no information about hepatic arterial pressure is available for those patients, pressure values at the entrance of hepatic artery were taken from the literature [START_REF] Ficher | The Liver[END_REF]. Pressure values at macro-cells level

Application to microspheres distribution in liver

Since our final aim is to simulate the injection of Theraspheres, we analyse their distribution under simplifying hypotheses, i.e. we simulate microspheres distribution as to be proportional to blood flow given by Poiseuille's law. Under such assumptions we are able to compute, depending on catheter tip placement, the percentage of injected healthy and tumour tissue, as well as the percentage of the injected dose which deposits in every tumour or healthy macro-cell. As shown in Fig. 4.21, we can propose an optimal injection point specific to every patient. In this context, optimisation refers to a cost function depending on injected healthy and tumour tissues percentages. The point suggested in Fig. 4.21a

(left) for Patient 1 induces radiations to 71% of the tumour and to 3% of the healthy liver. catheter tip is placed at the root of the arterial tree. In this case, all the tissue is irradiated, but a stronger concentration of radiations can be noticed in tumour cells, compared to normal ones. We quantified that in this scenario 9.2% of the injected activity deposits in the tumour for Patient 1 and 27.8% for Patient 2. This preponderance of microspheres in the tumour can be explained by two reasons: tumour macro-cells are more demanding in blood flow, and vascularisation is denser within the tumour. Even if only two injection points for patient are compared, the model allows the user to easily navigate into the vascular tree while calculating in real time the percentages of dose distribution mentioned above. This appears to be a promising approach in quantification of dose distribution in liver, with the aim of injection mode optimisation.

Validation of microspheres injection simulation

In order to estimate the validity of the simulation, we considered the SPECT-CT acquisition realised on Patient 2 after the 90 Y loaded microspheres injection. Table 4.5 shows the percentage of the injected activity that deposits to the tumour. The table compares the values given by the simulation to ground truth, calculated by an experienced nuclear radiologist of CEM using the "volumetric analysis" software (Syngo workstation; Siemens) and the method described in [START_REF] Garin | Boosted selective internai radiation therapy with 90 Y-loaded glass microspheres ( B-SIRT ) for hepatocellular carcinoma patients : a new personalized promising concept[END_REF]. In both cases and for both patients, the catheter was placed by the radiologist at the root of the tree. The discrepancy concerning the dose distribution for Patient 2 can be justified by multiple reasons. First, the patient presents several small tumour regions which were not segmented, and therefore not simulated, in this study. Second, the patient underwent a left lobe hepatectomy, which causes activity deposition in the tissues neighbouring the resection. Furthermore, clinical data concerning dose distribution were acquired through a software which performs tumours segmentation using a threshold applied to radiations intensity (cf. Fig. 1.7). This can lead to a biased computation of tumour volume, in particular to an overestimation of it, due to extra-tumour radiations deposition: a certain amount of healthy tissue surrounding the tumour can be considered as malignant if radiation intensity in such region is as high as in the tumour. In this case study, microspheres were injected at the root of the arterial tree. Blue-yellow colormap gives the percentage of microspheres in the injected tissue Fig. 4.22 shows a comparison between the SPECT-CT acquisition and the microspheres distribution simulation for the same patient 2. It can be observed that in both images microspheres concentration in tumour tissue is higher. This is due to the greater need of blood flow, characteristic of tumour tissue, and to the dense and tortuous architecture of tumour blood vessels. A precise comparison of liver shape, however, can not be achieved due to the nature of the two images. Indeed, SPECT-CT is composed of a stack of slices of the hepatic volume, whereas our model offers a 3D visualisation of the volume. The presented results were obtained starting from a simulated tree whose inlet diameter appears to be smaller than reality. Nevertheless, an input value of the model that will be unchanged in the improvement of such issue, is the ratio between blood flow in healthy and tumour tissues.

Discussion and conclusion

A hepatic arterial tree that mimics the contour and function of the patient's vasculature is generated. Results show good correspondence between real and in silico arterial trees.

This process requires segmentation of the shape of liver, tumours, and biggest and more visible hepatic arteries. The choice of the Frangi's parameters being relevant, we identified a set of parameters for each modality that appears to be suitable for any patients. In fact, the geometrical characteristics of the vessels (size, brightness) that Frangi's algorithm calculates strongly depend on the image modality, rather than to the single patient. For the same reasons we can deduce the stability of the parameters identified for the optimisation algorithm used for registration. In this context, a precise representation of vessels size is not necessary, due to the skeletonisation step that cancels this information. A more precise tailoring of Frangi's parameters will be needed when the actual vessel geometry will be used as a domain of fluid simulations.

Pressure value at the entrance of the hepatic artery, needed for the initialisation of the vascular model, is determined from the literature [START_REF] Ficher | The Liver[END_REF][START_REF] Lautt | Hepatic vasculature: a conceptual review[END_REF]. At macrocells level, pressure is estimated by general graphs concerning all the human systemic circuit (cf. Fig. 1.5). Based on such graphs, pressure values at macro-cells level can vary between 10 and 70 mmHg. Precise pressure values are difficult to obtain and very few studies are available in the literature [START_REF] Borowski | Pressure wire assessment of hemodynamic alterations after chemoembolization of hepatocellular carcinoma[END_REF]. It can be noticed that the current choice of pressure values provides a diameter of the entrance of the simulated hepatic artery that is smaller, for both patients, than reality. The pressure value at the macro-cells output may therefore be too small, or the pressure drop from the input to the output of the arterial tree may be underestimated in the model. Indeed, the pressure drop is currently considered only along each vessel segment, and not at the bifurcations relying two vessel segments.

The size difference can also be justified by a wrong estimation of the total amount of blood flow. Indeed, the optimisation criteria used during the vascular growth aim to minimise the total volume of the tree vessels. In reality, the vessels geometry may be more tortuous and induce a higher total volume than in the model. Moreover, both patients considered in this study present more than one small tumours that were not segmented and that may increase the total amount of blood flow. The input parameters of the model as pressure and blood flow will have to be adjusted in the future, based on patients' peculiarities. In particular, we plan to measure blood pressure at different levels of the hepatic arterial tree, thanks to a micro catheter. This subject will be discussed again in Chapter 5.

Some differences also appear between the geometry of the segmented arteries, used to initialise the model (Fig. segmented vessels is generally maintained. Nevertheless, using a rigid registration implies that some vessels irrigating peripheral regions of the liver may be poorly registered. This is the case for the vessel on the bottom left part of the tree in Fig. 4.23a. Indeed, the geometry of this vessel is maintained but, since it does not irrigate the liver, during the simulated vascular growth its radius progressively decreases. This is the reason why it is not visible in Fig. 4.23b, where, for the sake of clarity, only the vessels whose diameter is larger than 0.4 mm are shown.

Future improvements of the model will take into account the peculiarities of the vasculature due to different pathologies, like vessels tortuosity or arteriovenous shunt.

The algorithm defining the binary structure from segmented arteries will also need to be improved. Indeed, being the vessels skeleton very complex, it may happen that some of the tree branches are not detected. This can be due to a wrong treatment of bifurcation voxels clusters, or to the presence of trifurcations. The latter will be considered as a double bifurcation. Moreover, a drawback of this method is the lost of vessels diameter.

This should be integrated in the model, even if for the moment vessels radius can not be given as an initial parameter to the vascular model. Radius is indeed automatically defined by the model based on pressure and blood flow conditions.

We are aware that the simulated part of the arterial tree represents the real one only from a statistical point of view. Nevertheless, numerical simulation of microspheres transport and deposition still leads to a statistically realistic microspheres distribution map in the liver tissue, useful in clinical routine. A mean dose distribution map can be estimated by a kernel smoothing of the discrete distribution corresponding to the microspheres centers scatterplot. Moreover, such estimated dose distribution map will be nearly invariant to random changes of microspheres positions, as long as their statistical spatial distribution is preserved.

In order to better take into account liver functionalities, in the next chapters we will focus on the simulation of the blood flow and of microspheres trajectories simulation. It was noted [Childress and Kleinstreuer, 2014b] that microspheres transport could not be deduced directly from flow computation. The workflow we proposed here will be improved by recomputing flow and microspheres transport more accurately using CFD methods (cf. Fig. 1.9). Two kinds of simulation have been envisaged: a precise CFD simulation of particles trajectories down to the smallest vessels where a realistic mesh can be defined, and a simplified spheres distribution quantification as the one proposed in Section 4.3.3.

In order to couple the two models, the vascular growth model will have to be adapted in order to match with the CFD domain.

Simulation of the complete arterial tree has two more interests. First, pressure values given by Poiseuille's law at every level of the simulated tree can be used in order to initialise the pressure boundary conditions necessaries for the CFD simulations. Moreover, it can be used in order to perform a sensitivity study of blood flow, depending on boundary conditions on downstream vessels. Indeed, such conditions are known to be one of the principal difficulties in physiologically based CFD [START_REF] Olufsen | Modeling of the Arterial System with Reference to an Anesthesia Simulator[END_REF]. Complex and large vascular structures like those presented in this chapter will let us realise a sensitivity study of dose distribution, depending on catheter tip positioning.

Chapter 5

Patient-specific modelling of blood flow

When 90 Y microspheres are injected into the hepatic artery by the radiologist, the therapy takes advantage of the flow of blood, which transports the microspheres from the catheter injection position downstream in the systemic circuit, until they get stuck in arterioles or capillaries. The principal concern is to evaluate, as precisely as possible, how the microspheres are distributed in tumour and healthy tissues. Therefore, a sufficiently precise simulation of the flow field is needed before addressing the transport of microspheres in order to assist the radiologist.

Different factors can be at the origin of a poor treatment. It is well known that the presence of a tumour has an influence on the circulation of blood flow [START_REF] Oktar | Doppler sonographic evaluation of hemodynamic changes in colorectal liver metastases relative to liver size[END_REF][START_REF] Leen | Detection of hepatic metastases using duplex/color Doppler sonography[END_REF], related to the angiogenetic process [START_REF] Folkman | Tumor angiogenesis: therapeutic implications[END_REF][START_REF] Ledzewicz | Optimal and suboptimal protocols for a mathematical model for tumor anti-angiogenesis in combination with chemotherapy[END_REF].

Some tumours are strongly demanding in blood flow, others are less. In some particular circumstances, reflux can appear, meaning that blood flow, and thus microspheres, may move upstream in the systemic circuit. All these reasons justify the need for a detailed, patient-specific, 3D simulation of blood flow in the complex geometry of the hepatic artery, subjected to important inter-patient variability.

Clearly, different simulation approaches can be considered, from the more ambitious to very simple ones. A compromise between computational complexity and expected phenomena simulation accuracy must though be found. In Chapter 4 we presented a simplified approach: the principal arterial tree is extended by a simulated, complementary vessel network, which statistically reproduces the tumour and healthy tissue perfusion. Such simulation, based on Poiseuille's law, does not provide detailed fluid streamlines, which can actually be important for the particles transport. Indeed, the real streamlines can be tortuous in certain parts of the arterial tree. This leads to 3D simulations of the Navier-Stokes PDEs, on the domain delimited by the interior wall of the vessel tree, its inlet and its outlets. The choice of the boundary conditions (BC) on the domain borders represents a prerequisite that can be difficult to satisfy.

In this chapter, we will present and compare two of the principal methods displayed in the literature for the simulation of blood flow: the Finite Volume Method (FVM) and the Lattice Boltzmann Method (LBM). Some of the most common methods to circumvent the hindrance concerning BC and validation are reported, and a method based on MRI flow measures is proposed.

Another important hurdle consists in the impossibility to perform a full 3D simulation in a large vascular network, due to prohibitive computational costs. As shown in the next section, an alternative to full 3D simulation must be found in order to represent in a more numerically efficient manner all the sections of the systemic circuit where a precise description of the flow field is not needed.

Brief state of the art

Blood flow simulation is a very complex and wide research topic. Human blood consists in a suspension of cells in plasma, which is a Newtonian fluid mostly made of water. Blood mechanics is strongly influenced by the presence of red blood cells, especially in small vessels where their size is comparable to vessel size. Under such conditions, blood flow behaviour is non-Newtonian. Generally, blood flow is described by the Navier-Stokes equations, or a simplified version. Different models have been proposed in order to describe the non-Newtonian character of blood and its variable viscosity, among which the power law model and the widespread Quemada model [START_REF] Buchanan | Rheological effects on pulsatile hemodynamics in a stenosed tube[END_REF]. Other methods propose a precise description of blood cells as suspended particles in plasma [START_REF] Bessonov | Methods of Blood Flow Modelling[END_REF].

Blood pulsatility and the elastic response of arteries must also be considered for a more realistic simulation of blood flow when these phenomena must not be neglected.

As an alternative to the Navier-Stokes 3D approach, simpler numerical methods, developed for the solution of the hemodynamics equations, can be classified according to the domain dimensions. A review of zero and one dimensional models has been proposed

by [START_REF] Shi | Review of zero-D and 1-D models of blood flow in the cardiovascular system[END_REF]. Their advantage lies mainly in their strongly reduced computational cost, and the consequent possibility to enlarge the domain of the simulation even up to a simplified representation of the entire systemic circuit [START_REF] Bessonov | Methods of Blood Flow Modelling[END_REF]. As shown below, another noticeable advantage is that such methods can provide realistic downstream boundary conditions for 2D or 3D models to which they can be coupled.

Zero dimensional models, also called lumped parameter models, are based on an analogy between an hydraulic and an electrical circuit [START_REF] Smith | Physical Audio Signal Processing[END_REF][START_REF] Shi | Review of zero-D and 1-D models of blood flow in the cardiovascular system[END_REF][START_REF] Catanho | Model of Aortic Blood Flow Using the Windkessel Effect[END_REF][START_REF] Kokalari | Review on lumped parameter method for modeling the blood flow in systemic arteries[END_REF]. A certain section of the blood network is simplified as an electric circuit. In this context, resistance, capacitance and inductance of the electric circuit correspond, respectively, to frictional loss, vessel wall elasticity and blood inertia.

A model can be more or less precise based on the number of different compartments that are used to model the same network section. Typically, different compartments can be considered for the modelling of right heart, left heart, pulmonary circuit, and systemic circuit. The two latter can be further divided for more precision: aorta, arteries, arterioles, capillaries, and vice versa in the venous return. Each of these sections is known to have different properties concerning pressure, velocity, elasticity, pulsatility etc [START_REF] Martini | Fundamentals of Anatomy and Physiology[END_REF].

One of the most common lumped parameter model used in cardiovascular modelling is the Windkessel model, whose mathematical formulation was given by the german physiologist Otto Frank in 1899 [START_REF] Catanho | Model of Aortic Blood Flow Using the Windkessel Effect[END_REF]. Its simplest version considers all the venous network having zero pressure, and describes the arterial one, from aorta to capillaries, as a single resistor-capacitor circuit. Arterial compliance is described by the capacitor, and the resistance of the peripheral vessels by the resistor. In the three-elements Windkessel model, the resistance due to the aortic valve is also considered by a resistor, and in the 4-elements one an inductor, representing blood flow inertia, is added. The three or four elements Windkessel models are frequently used in multi-scale simulation of the cardiovascular circuit, for the definition of the downstream boundary conditions that can not be measured on patients, [START_REF] Vignon-Clementel | Outflow boundary conditions for 3D simulations of non-periodic blood flow and pressure fields in deformable arteries[END_REF][START_REF] Kung | In vitro validation of finite-element model of AAA hemodynamics incorporating realistic outlet boundary conditions[END_REF][START_REF] Bollache | Numerical modeling of arterial pulse wave propagation to characterize aortic hemodynamic: Validation using magnetic resonance data[END_REF][START_REF] Umbarkar | Computationally efficient fluid-particle dynamics simulations of arterial systems[END_REF], Childress and Kleinstreuer, 2014b] among others.

One-dimensional models well describe the pressure wave propagation in arterial circuit, but still provide averaged information on the blood pressure and flow field, thus not precise enough if the simulation of transported microspheres trajectories is to be taken into account. Nevertheless, the obtained averaged fluid results can be used in more complex 3D models, for the definition of more realistic boundary conditions distally in the arterial tree, where in vivo flow values are impossible to measure. [START_REF] Formaggia | One dimensional models for blood flow in arteries[END_REF] propose a one dimensional model based on a system of partial differential equations modelling the evolution of pressure wave propagation in compliant arteries. The Navier-Stokes equation is integrated over the flow cross-section and is coupled to a mechanics equation describing vessels wall displacement. [START_REF] Olufsen | Numerical simulation and experimental validation of blood flow in arteries with structured-tree outflow conditions[END_REF] propose a one dimensional description of blood flow in all the principal systemic arteries, where the peripheral arterial network is considered as a structured, asymmetric binary tree. The structured tree model for the simulation of peripheral arteries provides physiological boundary conditions to a model of blood flow in large arteries based on a linearised Navier-Stokes equation. In [START_REF] Bollache | Numerical modeling of arterial pulse wave propagation to characterize aortic hemodynamic: Validation using magnetic resonance data[END_REF] a fluid-structure model of blood flow in aorta is proposed. The 1D model is coupled to a Windkessel model of peripheral arteries and validated against phase contrast MRI.

Multidimensional models of blood flow. In case of axisymmetric tubes, a two dimensional simulations can provide a description of the radial velocity profile. Nevertheless, full 3D models are often needed for more complex domains and velocity profiles. Three dimensional computational fluid dynamics allows accurate simulation of blood flow field and its interactions with vessel walls (fluid-structure interaction models). The literature proposes numerous studies on CFD simulations of blood flow in the cardiovascular system.

In particular, since the mid-1990s, researchers started focusing on patient-specific simulations, due to the possible integration of medical imaging in CFD [START_REF] Steinman | Flow imaging and computing: Large artery hemodynamics[END_REF]. The finite elements [START_REF] Taylor | Finite element modeling of blood flow in arteries[END_REF]] and finite volumes methods are often used for the solution of the incompressible Navier-Stokes equation and 3D models are often coupled to a zero or one dimensional model for the simulation of input and output impedances [START_REF] Vignon-Clementel | Outflow boundary conditions for 3D simulations of non-periodic blood flow and pressure fields in deformable arteries[END_REF].

Generally, in vivo measurements of blood velocity are performed based on phase contrast MRI or Doppler US and blood pressure can also be measured through specific devices,

as in [START_REF] Kung | In vitro validation of finite-element model of AAA hemodynamics incorporating realistic outlet boundary conditions[END_REF]. Such values are used for initialisation and/or validation of numerical simulations. In [START_REF] Kung | In vitro validation of finite-element model of AAA hemodynamics incorporating realistic outlet boundary conditions[END_REF] Multidimensional models of hepatic arterial blood flow. A similar procedure applied to the hepatic blood flow has been less investigated. The smaller vessels size and the complexity of the geometry make it harder to extract informations concerning vessels geometry and flow field. As already mentioned in Chapter 1, patient-specific simulations of blood flow in hepatic arteries have been proposed by [START_REF] Xu | Solid Tumor Embolotherapy in Hepatic Arteries with an Anti-reflux Catheter System[END_REF] and Aramburu et al. [2016a]. The former research team proposes a simulation domain consisting on a 3D but planar representative hepatic arterial system with 17 outlets. The finite volume method is used here, through the Ansys-Fluent software, for 3D particle-hemodynamics simulations, where blood is modelled as a laminar non-Newtonian fluid following a modified Quemada model [START_REF] Buchanan | Rheological effects on pulsatile hemodynamics in a stenosed tube[END_REF]. Our contribution, presented in Sections 5.2.2 and 5.3.3 of this chapter, proposes a comparison between the classical finite volumes method and the lattice Boltzmann method.

The computational domain of the presented simulations performed in Ansys Fluent is the actual patient's geometry, where no simplifications about the shape of vessels section is made, nor on the symmetry of the bifurcations. The method based on LBM is developed in order to be able to simulate blood flow in the hepatic arterial tree generated with the method described in the previous chapter.

Lattice Boltzmann method

The Lattice Boltzmann method derives from the statistical kinetic theory of gases, where the Boltzmann equation describes the evolution in time, based on micro-dynamic interactions, of the non-equilibrium (i.e. depending on time) distribution density f (x, p, t). The The Boltzmann equation is an integro-differential equation, whose discretised version, the lattice Boltzmann equation, is numerically solved on a regular lattice in R 3 [Succi, 2001, Sukop and[START_REF] Sukop | Lattice Boltzmann Modeling: An Introduction for Geoscientists and Engineers[END_REF]. In case of incompressible fluids with low Mach number, it has been shown that the Navier-Stokes equations can be deduced at the second-order of approximation from the lattice Boltzmann equation through a Chapman-Enskog expansion [START_REF] Qian | Lattice BGK Models for Navier-Stokes Equation[END_REF], He and Luo, 1997b[START_REF] Chen | Lattice Boltzmann method for fluid flows[END_REF][START_REF] Luo | Lattice Boltzmann Methods for Computational Fluid Dynamics[END_REF][START_REF] Wang | Evaluation of three lattice Boltzmann models for particulate flows[END_REF]. We recall that the Mach number is a non-dimensional quantity defined as Ma = v cs , where v is the characteristic velocity of the fluid and c s the speed of sound in the fluid.

Such value describes the balance between the inertial effects and the pressure gradients. It can be noted that Boltzmann equation and lattice Boltzmann equation can be also defined in 2D, for (x, p, t) ∈ R 2 × R 2 × R.

In the classical lattice Boltzmann method, the simulation domain is discretised in a regular lattice characterised by a fixed distance h between two adjacent nodes along the x and y and, in the 3D case, z directions.

The fluid to be simulated can be considered to be made of virtual, mesoscopic particles, whose positions are constrained to belong to the set of lattice nodes. At every time step, each particle situated on a certain node is constrained to move along a straight line towards one its neighbours. The set of authorised neighbouring nodes includes the initial one and depends on the discretisation scheme.

Hence, in the lattice Boltzmann context, the phase space is discretised and the cubelet ∆x is considered as the cube centred on a certain lattice node with a h linear size, and particles velocities e i are constrained to a limited set I of velocities, that depends on the chosen discretisation model. Therefore, e i ∈ I are the authorised elementary displacement.

It can be imagined that particles stream from a node to another, colliding with each other according to predefined deterministic rules prescribed to respect mass, momentum and energy conservation. The algorithm does not move particles explicitly. Instead, the explicitly used value is a discrete equivalent of the distribution density f (x, p, t) appearing in the Boltzmann equation, namely the discrete distribution function f i (x, t).

This distribution gives, at the discrete time t, the expected number of particles starting from the node x on the lattice, and moving with the velocity e i ∈ I, during one time step of the LBM algorithm. The quantities f i (x, t) are updated applying simple rules, in an identical manner at every node x of the lattice. This is the reason leading to very efficient implementations on generic parallel processing architectures, as the graphical parallel unit (GPU) [START_REF] Jurczuk | GPU-based computational modeling of magnetic resonance imaging of vascular structures[END_REF].

The macroscopic quantities of the fluid are reconstructed at time t and at each node x as weighted averages of the mesoscopic distributions f i (x, t), i ∈ I. Thus, the fluid mass density and velocity at node x and time t are the zeroth and first moments of the probability density f i (x, t), i ∈ I, defined as:

ρ = i∈I f i , v = 1 ρ i∈I f i e i .
(5.1)

The LBM algorithm modifying the distributions f i (x, t) at every time step is considered below. In an ideal case where no collision occurs, the evolution of the distribution function must verify:

f i (x + e i ∆t, t + ∆t) = f i (x, t), (5.2)
where every particle residing on the node x at time t and having a velocity e i propagates to the node x + e i ∆t at the time t + ∆t, conserving the velocity e i . Due to the collisions occurring between particles, a collision operator must be considered in order to take into account all particles that arrive to the node x + e i ∆t with velocity e i but that at time t were elsewhere than the node x and viceversa, particles residing at node x at time t with velocity e i that at time t + ∆t do not reach the node x + e i ∆t.

Different models for the collision operator exist. We will present the most common one introduced by Bhatnagar Gross and Krook (BGK) [Bhatnagar et al., 1954, Chen and[START_REF] Chen | Lattice Boltzmann method for fluid flows[END_REF]] but others exist, like the Multiple Relaxation Time model [START_REF] D'humières | Multiple-Relaxation-Time Lattice Boltzmann Models in Three Dimensions[END_REF], the regularised model [START_REF] Latt | Lattice Boltzmann method with regularized pre-collision distribution functions[END_REF]] and the entropic model [START_REF] Ansumali | Minimal kinetic modeling of hydrodynamics[END_REF]. In the context of the single-relaxation-time BGK formulation, the evolution over time of the distribution function is defined as:

f i (x + e i ∆t, t + ∆t) = f i (x, t) - ∆t τ f i (x, t) -f eq i (x) i ∈ I, (5.3)
where τ is the relaxation time and f eq i is the local equilibrium distribution function which describes the equilibrium state when the particles gain and loss of a same node are perfectly balanced. The τ value is a dimensionless parameter that gives the rate of approaching to equilibrium. The equilibrium distribution functions are defined as (5.4) where C is the lattice speed defined as C = h ∆t . The weight factors w i derive from a low Mach number expansion of the Boltzmann-Maxwellian equilibrium distribution function and from the conservation laws. This procedure is not straightforward and details can be found in [He andLuo, 1997b, Succi, 2001]. The values of the weights w i in the most common 2D and 3D models will be given below.

f eq i = w i ρ 1 + 3 C 2 e i • v + 9 2C 4 (e i • v) 2 - 3 2C 2 v • v ,
In the lattice Boltzmann formulation, fluid variables like density and viscosity are not defined in an explicit way but are introduced indirectly in the algorithm. The relaxation time is strongly related to the fluid kinematic viscosity ν, which is computed as:

ν = 1 3 τ - 1 2 h 2 ∆t = 1 3 τ - 1 2
Ch.

(5.5)

It can be seen that in order to recover a physically realistic viscosity value, the relaxation time must verify τ > 1 2 . It has also been shown that a value close to 1 guarantees more computational stability [START_REF] Sukop | Lattice Boltzmann Modeling: An Introduction for Geoscientists and Engineers[END_REF]. On the other side, fluid density is related to pressure field p by the state equation

p = C 2 s ρ , (5.6)
where C s is the lattice speed of sound and can be deduced by the lattice speed through the equation

C s = C √ 3 .
(5.7)

A time step of the lattice Boltzmann algorithm is divided in two principal substeps:

streaming and collision. Equation (5.3) could actually be solved in one step. Nevertheless, due to the particular numerical formulation of boundary conditions on solid boundaries, it is a common, more practical approach to separate the two phases. In a first step, the particles velocities are updated according to the collision operator. Defining t + as the post-collision time, the following relationship is thus computed:

f i (x, t + ) = f i (x, t) - ∆t τ f i (x, t) -f eq i (x) i ∈ I, (5.8)
Then, during the streaming phase, the virtual particles (considered implicitly) either stay on the same node or jump to a neighbouring one according to their new velocity:

f i (x + e i ∆t, t + ∆t) = f i (x, t + ) i ∈ I.
(5.9) from [START_REF] Ehrhardt | Novel Trends in Lattice-Boltzmann Methods[END_REF] Discretisation of physical domain and velocity space Different models for the discretisation of the phase space exist, generally labeled as density and velocity of Eq. (5.1). Whereas the different components of the velocity vector can explicitly be defined in terms of the distribution functions, pressure is specified through the fluid density, thanks to the equation of state (5.6). Details of the implementation of the boundary conditions in the lattice formulation can be found in [START_REF] Zou | On pressure and velocity boundary conditions for the lattice Boltzmann BGK model[END_REF]He, 1997, Sukop and[START_REF] Sukop | Lattice Boltzmann Modeling: An Introduction for Geoscientists and Engineers[END_REF].

Despite its easy implementation and computational efficiency, such method is limited to computational domains whose boundaries are parallel to the lattice axes. A method able to handle curved boundaries was proposed by [START_REF] Guo | An extrapolation method for boundary conditions in lattice Boltzmann method[END_REF]. In this method the distribution functions f i of the boundary nodes whose velocities are directed towards the interior of the domain, are split in their equilibrium and non-equilibrium parts. An extrapolation method is thus used to deduce such distribution functions f i on the boundaries.

A particular boundary condition that is worth to describe in details (the reason for which will be clarified in Chapter 6) is the so called "bounce-back" rule, which enforces no-slip conditions on the walls of the domain. The method is based on the idea that the macroscopic fluid velocity is zero at the lattice boundaries if all mesoscopic particles approaching the boundaries are bounced back with the same but inverse velocity. More precisely, if the direction i links a node inside the domain to one that is outside, the distribution function of a particle crossing this link will be flipped and, at the following step, it will have the opposite direction -i:

f -i (x, t + 1) = f i (x, t + ), (5.14)
where t + denotes the time just after the collision but before the following streaming step.

Such rule satisfies the condition of zero tangential velocity along the computational domain.

In this method the boundary between the fluid and the wall is considered to be exactly halfway between the solid and the fluid nodes connected by the link i.

Conversion to lattice units

The numerical efficiency of the Lattice Boltzmann method is also related to the simplicity of the implemented equations. In particular, a different, non-dimensionalised, system of coordinates is used, which is based on the unity of temporal and spatial discretisation parameters: the lattice units being marked with an overbar, we define ∆t = 1, h = 1 and consequently C = 1. This procedure, though, complicates the recovering of physical values and, most importantly, the physical intuition.

Different procedures for the parameters selection have been suggested in the literature.

For example, [START_REF] Feng | Coupled lattice Boltzmann method and discrete element modelling of particle transport in turbulent fluid flows : Computational issues[END_REF] propose to define first h based on computational cost and τ in order to limit numerical instabilities. The simulation would thus be launched with the ∆t computed from (5.5) and, if convergence is reached, the Mach number can be checked to evaluate the accuracy of the solution found.

The coherence between the physical and the lattice units can be ensured through the Reynolds or the Mach number, which are both non-dimensional parameters. This can also help in the choice of the lattice parameters. In our C++ in-house LBM software we decided to follow a procedure similar to the one proposed by [START_REF] Feng | Coupled lattice Boltzmann method and discrete element modelling of particle transport in turbulent fluid flows : Computational issues[END_REF]. The parameter h is manually chosen based on the given domain, but τ is computed based on the Ma 1 condition. Based on our observations, a Mach number smaller than 0.05 is recommended.

Since this quantity is non-dimensional, the same value must be recovered if it is computed starting from physical or lattice units:

Ma = v max C = v max C = v max .
(5.20)

The conversion between lattice and physical velocities (5.18) gives:

C = v max v max .
(5.21)

If the maximum physical velocity is available (for example when considering Poiseuille flow), then the conversion parameter C can be computed in this way. Otherwise, an estimation should be found based on the domain geometry and boundary conditions. Once the spatial discretisation h and the lattice velocity C are known, the temporal discretisation ∆t can be computed. Finally, (5.5) gives a relaxation time τ value that is tailored on the considered simulation problem.

LBM blood flow simulation in a patient-specific hepatic artery

The method described above has been applied for the simulation of blood flow on the vessels geometry extracted from the patients data, according to the procedure described in Section 4.2. The entire simulated tree was cut in order to simulate blood flow only in the principal vessels, corresponding to the segmented ones. This let us reduce, in a first step, computational time, and compare the results to those obtained using the finite volumes method in the same part of the hepatic artery. Different diameter values were used to filter the simulated tree in order to select the largest vessels; some examples are given in Fig. 5.2. The initial conditions of the vascular model that are used in Chapter 4 provide a diameter of the entrance of the simulated hepatic artery that is smaller than reality. Such initial conditions were then adjusted in order to obtain the right diameter, at least at the inlet of the arterial tree. For both patients, a pressure value at macro-cells level of 80 mmHg led to a simulated arterial tree whose input diameter corresponds to the right patient's one. Under such conditions, the total pressure difference between input and macro-cells (90 mmHg and 80 mmHg) may seem very small, compared to physiological data. Actually, even if pressure value at capillaries is known to vary between, 2 and 10 mmHg, it is also known that 98% of the pressure drop occurs just before the capillaries [START_REF] Lautt | Hepatic vasculature: a conceptual review[END_REF]. Moreover, Borowski et al.

[2013] measured blood pressure, with a catheter placed distal in the hepatic arterial tree, as 71±26 mmHg. Such information, coupled to the very rough general systemic pressure graph (Fig. 1.5), shows the difficulties in identifying precise values. The obtained tree was finally cut in order to remove all the vessels with a diameter smaller than 2 mm. The resulting geometry was thus consistent with the segmented one.

The LBM simulation was then initialised using, as boundary conditions, the pressure values provided by the vascular growth model at every bifurcation of the simulated tree.

Concerning the blood simulation, important simplifications were applied to the arterial flow. In this first phase, we neglected blood pulsatility, as well as vessels compliance and the non-Newtonian character of blood in small arteries. Thus, blood was considered as an incompressible Newtonian fluid with density ρ=1051 kg/m 3 and constant kinematic viscosity ν=3.8 10 -6 m 2 /s [START_REF] Levick | An Introduction to Cardiovascular Physiology[END_REF].

The procedure described in [START_REF] Jurczuk | In Silico Modeling of Magnetic Resonance Flow Imaging in Complex Vascular Networks[END_REF] is used in order to convert the vessel segments of the vascular model into three dimensional vessels whose shape is described by truncated cones.

Results

The described procedure was applied to both patients already presented in Chapter 4. The LBM simulation was launched using, as the computational domain, the principal vessels of the geometry simulated using the vascular growth model. In this work, the LBM program was launched in its non-parallelised version, and the computation lasted between 3 and 5 hours. The used numerical parameters are: ∆t = 3.18 10 -6 s, h = 1.3 10 -4 m, τ = 0.94.

The LBM simulations of blood flow under the described characteristics for both patients are shown in Fig. 5.3 and Fig. 5.4. Due to no-slip conditions, blood flow velocity against the vessel walls is zero. The vessels are thus divided in regular slices for visualisation purposes.

Finite Volumes method

As most of the methods used for the discretisation of partial differential equations, the finite volumes method is based on the discretisation of the spatial domain in a computational mesh. In the FVM the domain is divided in the so called control volumes. The differential form of the governing equations is then integrated over each control volume. This method is particularly adapted to the solution of conservation equations like the Navier-Stokes equations, because it ensures the conservation of quantities like mass, momentum and energy, over each control volume. A general finite volume discretisation scheme for conservation laws and in particular for the Navier-Stokes equation is described in details in Appendix B.

Generation of a CFD mesh from a clinical image

The software Ansys Fluent, whose mathematical background functioning is described in Appendix B, was used for the FVM simulation of blood flow in a patient-specific geometry, The irregular segmented binary volume was converted into a regular meshed geometry, from which stable CFD simulations can be performed. The procedure is not straightforward and is addressed here.

In order to implement blood flow simulation with a CFD software like Ansys Fluent, a closed smooth mesh is needed. In such mesh, every different surface corresponding to the input, and to all the outputs of blood flow, must be identified as distinguished structures and independently imported in Fluent. The binary volume resulting from the Frangi-based segmentation is thus used to create a triangular mesh corresponding to the surface of the segmented vessels. This preliminary mesh is generated using the "contour" function in the Paraview software. Input and output surfaces selection is achieved by a manual cut of the mesh at every vessel extremity. The result is a mesh describing a surface which is open at the entrance of blood and at every vessel end. A number of surfaces, called cap, corresponding to the number of mesh holes, needs to be created in order to perfectly close the open mesh. Such caps can be created thanks to the properties of the mesh: indeed, mesh consists of only triangles. When the mesh is clipped by a plane, new mesh nodes are created on the intersection between the mesh triangles and the plane, transforming some of the triangles into quadrilaterals. We developed an algorithm that first looks for all the mesh nodes which belong to a triangle or a quadrilateral, and which belong to only three facets or less. Such property guarantees that they are on the cut border of the mesh. All the identified border points then need to be distinguished between every cap. A k-means clustering algorithm [START_REF] Arthur | K-means++: The Advantages of Careful Seeding[END_REF] accomplishes this step minimising the average squared Euclidean distance into every group. Cap surfaces are created by a simple preliminary triangularisation: for every set of nodes the barycentre is calculated, points are then classified in a trigonometric order, and triangles are finally defined by linking every point to its following neighbour and to the barycentre. An example of such procedure can be seen in Fig. 5.5. Finally, quadrilaterals resulting from mesh clipping are split in two triangles in order to redefine a proper triangular mesh. Afterwards, the surface corresponding to vessels wall and those corresponding to caps are imported in ANSYS Fluent meshing. A remeshing is performed using wrapping in order to obtain a smooth closed surface and mesh quality compatible with fluid simulations. Triangles minimum and maximum sizes can be adapted to the segmented vessels in order to smooth the surface or to better describe bifurcations and narrow vessels by choosing, respectively, a larger or smaller triangle size.

Different parameters for minimum and maximum sizes of the mesh elements were tested.

A minimum size of 0.3 mm and a maximum one of 2.5 mm were chosen and give a good, smooth, representation of the geometry. In order to capture the behaviour of the boundary layer with sufficient precision and reduce numerical instabilities, multiple layers of prisms are capping the artery walls, which is a common approach when using a finite volume solver as Fluent. A volume mesh is automatically generated on the basis of the surface mesh and passed to the ANSYS Fluent Solver described above.

FVM blood flow simulation in a patient-specific hepatic artery

In this context as well, blood is preliminarily considered as an incompressible Newtonian fluid. Its simulated dynamic viscosity value is defined as µ = 0.004 Pa s, and its density as ρ = 1060 kg /m 3 [START_REF] Levick | An Introduction to Cardiovascular Physiology[END_REF]. Inlet condition at the proper hepatic artery is given at mean arterial blood velocity, which is defined to be 30 cm/s over the entire vessel section based on literature values [START_REF] Dyvorne | Quantification of hepatic blood flow using a highresolution phase-contrast MRI sequence with compressed sensing acceleration[END_REF][START_REF] Tziafalia | Echo-Doppler measurements of portal vein and hepatic artery in asymptomatic patients with hepatitis B virus and healthy adults[END_REF][START_REF] Hubner | Hepatic arterial blood flow velocities: assessment by transcutaneous and intravascular Doppler sonography[END_REF][START_REF] Kito | Doppler sonography of hepatic arterial blood flow velocity after percutaneous transhepatic portal vein embolization[END_REF]. In such studies blood flow is measured in either the common, left or right hepatic artery, and the average velocity varies between 11 and 50 cm/s. Output boundary conditions are given on blood flow percentage, as to be proportional to every outlet vessel diameter. Classical no-slip conditions are enforced on vessel walls.

Steady-state Navier-Stokes equations:

∇ • u = 0 (5.22) u • ∇u = - 1 ρ ∇p + ν∆u (5.23)
are solved in ANSYS Fluent. A pressure-based solver was used together with a segregated SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) scheme (cf. Appendix B).

The convergence is verified by checking the mass and momentum conservation in between inlet and outlets. Calculation are stopped when a difference of 10 -7 for continuity and 10 -4 for momentum is reached. To enhance convergence, a Least-Square Gradient Method was chosen. Such method is comparable to the one of the node-based gradient on skewed and distorted unstructured meshes and is less expensive to compute [START_REF] Rausch | Spatial Adaptation of Unstructured Meshes for Unsteady Aerodynamic Flow Computations[END_REF].

Results

The procedure was applied to the cone-beam angiographies of the 2 patients considered before plus 2 other patients suffering from an advanced stage HCC, as can be seen in Moreover, metallic clips near the patient's tumour make the image particularly noisy and complicate the segmentation of the tumour vessels. Vessels segmentation results in arterial trees with 20 to 50 outlet vessels. Radius of the proper hepatic artery varies between 2.3 and 3.6 mm. The smallest radius of segmented vessels is 0.49 mm. Details about the obtained mesh are given in Table 5.1 for all the patients. Min. face area (mm 2 ) 2.5 e -3 1.1 e -3 1.9 e -3 1.1 e -3

Max. face area (mm 2 ) 0.56 0.298 0.589 0.47

Validation

The problem of the validation of blood flow simulation is an open issue in the literature.

The implementation of realistic simulations mostly relies on the choice of the right bound-ary conditions for the initialisation of the simulation, as well as on the definition of a domain which well represents the vessels shape. Indeed, commercial softwares like Ansys

Fluent have extensively been validated by the vendor, and our in-house LBM software has already been validated as shown in previous papers [START_REF] Jurczuk | Computational modeling of MR flow imaging by the lattice Boltzmann method and Bloch equation[END_REF][START_REF] Jurczuk | In Silico Modeling of Magnetic Resonance Flow Imaging in Complex Vascular Networks[END_REF]. We can therefore state that, if the right boundary conditions are chosen, the simulations results will consequently be correct, depending on the convergence error.

A comparison of the results obtained using the two different methods was also realised.

The comparison was possible due to the similar geometries. Indeed, as already mentioned in Section 5.2.1, the arterial tree simulated with the vascular growth model was initialised with pressure values providing a correct diameter of the tree inlet. For flow simulation, the choice of the boundary conditions was different: velocity was used in the FVM simulation, and pressure for the LBM simulations. Nevertheless, the pressure values provided by the Poiseuille's law at every bifurcation of the simulated tree were consistent with the values computed by the Fluent simulations. The LBM simulations were thus initialised using pressure conditions very similar to the values computed by the FVM. Indeed, the pressure differences between inlet and outlets varied, in both cases, approximately between 1 and 5 mmHg. It can be seen in Fig. 5.10 and Fig. 5.11 that the obtained velocities are qualitatively comparable. Some differences naturally appear due to the differences in the geometries of the computational domains, deriving from different procedures. Inaccuracies in the LBM simulation may derive from the simplification of the vessels geometry to a set of regular, truncated cones, that do not correspond exactly to the extracted vessels because their diameter is modified during the vascular growth. Even if the centrelines of the extracted vessels are unchanged, the vessels diameters are adjusted based on the amount of liver volume that they irrigate and on the output pressure. Therefore, when only the vessels whose diameter is larger than 1 mm are selected for the simulation, some of the input vessels may be excluded. Possible vessels spasms are also neglected in the LBM simulation: the vessels diameter is not an initialisation parameter, but is automatically recomputed by the vascular growth model. The LBM simulation is also influenced by the ratio between blood flow in healthy and in tumour macrocells. This value is usually considered to be at least 2 but can be much higher [START_REF] Andrews | Modulation of liver tumor blood flow with hepatic arterial epinephrine: a SPECT study[END_REF]] and strongly influence the blood behaviour. The influence of such parameter, defined in this work as 3, on the geometry of the simulated tree and on the obtained velocities will have to be explored. On the other side, the FVM simulations in the patient-specific arteries are subject to vessels segmentation in-accuracies. In particular, the diameter of the smaller vessels tends to be overestimated due to the non adaptive thresholding of the image after the Frangi's enhancement. This may slow down the blood simulated velocity. Vessels irrigating the tumours are also difficult to distinguish from one another due to their tortuosity, and may be erroneously segmented as one, larger vessel. In this work, an image-based approach is used to estimate as much information as possible from the patients clinical images. As already shown in Chapter 2, we performed 3D PC MRI phase contrast acquisitions on healthy volunteers and HCC patients, which gave the results shown in Table 2.5 Sec. 2.2.2, concerning mean blood velocity until a vessel diameter of approximately 2 mm. The pixel size was here defined as 0.9 mm.

Due to clinical constraints, it was not possible to perform PC MRI on the same patients as those whose hepatic arteries were segmented. A representative value of velocity was then used in these simulations and a precise validation was therefore impossible. Nevertheless, we observe a coherence between the results obtained from the simulations and the velocity measurements realised on PC MRI.

In order to support the choice of the output boundary conditions of the simulations, the ratios between diameter and blood flow in vessels of similar sizes were studied on a 3D

PC MRI sequence acquired on a healthy volunteer (Table 5.3). It is well known that vessels PC MRI were possible only on vessels much larger than the outputs of the simulation, we verified that our choice of boundary conditions was reasonable. The ratios between diameter and blood flow are similar in all considered vessels, except in the left hepatic artery.

Nevertheless, the diameter of this vessel is smaller and, consequently, the measurement less reliable.

The results obtained with the simulations concerning pressure and velocity in larger arteries are coherent with literature [START_REF] Ficher | The Liver[END_REF] and with our measurements (cf. Table 2.5). Further in the hepatic arterial tree, approximate information relying pressure and velocity to vessels diameter in human systemic circuit exists. Pressure of arteries at sizes as those analysed in current study can vary approximately between 70 and 100 mmHg [START_REF] Martini | Fundamentals of Anatomy and Physiology[END_REF][START_REF] Borowski | Pressure wire assessment of hemodynamic alterations after chemoembolization of hepatocellular carcinoma[END_REF][START_REF] Bura-Rivière | Physiologie et exploration de la circulation artérielle[END_REF]]. Information about blood velocity at small hepatic arteries is also hard to determine exactly. According to [Bura-Rivière and Boccalon, 2010], mean blood velocity for arteries whose diameter varies between 2 and 6 mm, can vary between 20 and 50 cm/s.

Conclusion

We developed a procedure to extract a patient-specific hepatic arterial tree from cone-beam CT 3D angiography, until a vessel diameter of around 0.5 mm. We realised blood flow simulation using the lattice Boltzmann and the finite volumes methods. The simulation results concerning blood pressure and velocity appear to be coherent with the literature, even if flow precise characteristics at such level are very little known and vary significantly between different individuals and under different circumstances.

In this study, the LBM simulation was launched in a non-parallelised version and the computational time was noticeably longer than the parallelised FVM simulation, on a similar geometry. Nevertheless, the LBM flow computation in the simulated arterial tree using a computer cluster was shown to be between 1.7 (2 processors) and 47 (128 processors)

times faster [START_REF] Jurczuk | In Silico Modeling of Magnetic Resonance Flow Imaging in Complex Vascular Networks[END_REF]. Using a single GPU, speeds up the computation of 10-100 times, as shown in [START_REF] Jurczuk | GPU-based computational modeling of magnetic resonance imaging of vascular structures[END_REF]. Based on the experience in the research time in parallelisation, an implementation of the code for multi-GPU is planned and we can expect to reduce the computational time by a factor of some hundreds. This will let us simulate the blood flow in a much more complex vascular structure, as the one presented in Fig. 5.12 or even more complex. Our MRI measurements are not statistically representative but are validated by an experienced radiologist and by known mean velocities at the aorta and inferior vena cava levels [START_REF] Gabe | Measurement of Instantaneous Blood Flow Velocity and Pressure in Conscious Man with a Catheter-Tip Velocity Probe[END_REF]. Their interest is double: information about blood velocity at smaller arteries is crucial for the definition of boundary conditions, by letting us tuning output conditions that are currently little known, in particular tumours blood flow. Output conditions should be adjusted in order to obey to some data fidelity relative to velocity measured on PC MRI. Moreover, PC MRI will let us validate blood flow results when it will be available on the same patients as 3D angiography.

New 4D PC MRI techniques

, not yet available on the MR system used for this study, should allow time-resolved acquisitions of three dimensional PC MRI, synchronised to the cardiac cycle [START_REF] Stankovic | 4D flow imaging with MRI[END_REF]. It should be noticed though, that time constraints in the clinical workflow do not actually encourage this approach.

A more precise description of variable blood viscosity and pulsatility will also be performed. Moreover, before SIRT, patients can undertake chemoembolisation, vessels coiling or clipping whose metal can induce artefacts in acquisitions and reduce the quality of the segmentation, as for Patient 1, where the diameters of tumour vessels near the clips are overestimated. Numerical results are inevitably dependent on segmentation, which 178CHAPTER 6. PATIENT-SPECIFIC MODELLING OF MICROSPHERES TRANSPORT this class can be defined by a particle size being smaller than 1 to 10 µm. In this case, the liquid with particles is called a colloidal dispersion, and can be numerically addressed as a continuum. Simulations of fluidised beds, for example, can be fully resolved by methods where particles are considered to be a solid phase of a non-Newtonian, compressible fluid (e.g. the Two-Fluids Model, [START_REF] Kuipers | A numerical model of gas-fluidized beds[END_REF]).

The intermediate size class is characterised by the dominance of viscous forces on Brownian forces. As shown by [START_REF] Stickel | Fluid Mechanics and Rheology of Dense Suspensions[END_REF], suspensions falling into this category can be simulated by continuum methods. However, this might not give a satisfactory representation of the rheological behaviour [START_REF] Leonardi | Coupled DEM-LBM method for the free-surface simulation of heterogeneous suspensions[END_REF], and a more precise representation of particles might be preferred.

The last class is made by larger particles, where the collisional forces between them are the dominant forces.

Another important feature should be analysed in order to understand the nature of the two-phase fluid: the density of solid grains present in the liquid. Many definitions of a "dense" suspension exist. A suspension can be considered as "highly concentrated" if the solid volume fraction is greater than 50% [START_REF] Ladd | Lattice-Boltzmann Simulations of Particle-Fluid Suspensions[END_REF]]. We will consider a suspension "dense" when the average distance between grains is smaller than the average grain size [START_REF] Stickel | Fluid Mechanics and Rheology of Dense Suspensions[END_REF]. A "dilute" suspension is characterised by a negligible solid fraction φ, which is defined in [START_REF] Succi | The Lattice Boltzmann Equation for Fluid Dynamics and Beyond[END_REF] by: φ = 4π p a 3 p 3L 3 , (6.1) a p being the radius of the particle p and L the size of a 3D box containing the fluid. More precisely, Ladd [1994b] shows that three typical conditions exist: a suspension is considered to be concentrated if φ is close to 0.45, semi-dilute if φ 0.25, and dilute if φ 0.05.

The simulation of the transport of solid microspheres can be performed through many numerical methods, where the particles can be considered as a continuum or not. In this work, we will concentrate on two principal families of methods: the Lattice Boltzmann Method and the traditional Computational Fluid Dynamics (in particular, the Finite Volumes Method).

the classical bounce back rule. In that case, the bounce-back rule is used to translate into the Lattice Boltzmann context the classical no-slip condition on the stationary boundaries of a computation domain (cf. Chapter 5). In this context, the phenomenon is very similar except that the boundaries (here, the particles) are moving with the fluid. The criterion used in this context is indeed called the "modified bounce-back rule". Equation (5.14) is hence revised in order to take into account the velocity of the particle.

Being the diameter of Theraspheres ® small enough (∼ 20-30 µm), we decided to neglect the rotation of the spheres, and to consider therefore that the velocity of the particle is constant on all its nodes.

The modified bounce-back rule is defined as:

f -i (x, t + 1) = f i (x, t + ) -6w i ρv p • e i , (6.2)
where i indicates a boundary link and -i the opposite one, as shown in Fig. 6.2 for the two dimensional case. As already mentioned in Chapter 5, f i (x, t + 1) is the post-collision density distribution at the node x, whereas w i are the same weighting factors as in the equilibrium distribution function, and ρ is the fluid density at the node x. The variable v p represents the velocity of the particle p, and e i the ith LBM elementary direction. As shown in Fig. 6.2, equation (6.2) means that the LBM mesoscopic particles that approach the boundary links are reflected back with their same velocity, minus a value proportional to the velocity of the particle.

A peculiarity of this method is the fact that both fluid and solid boundary nodes are considered in the same way. Indeed, equation (6.2) is applied to the boundary nodes of both natures. This method indeed considers the particle as a sort of "shell", filled with fluid. This clearly does not respect the physics of a solid particle, but it has been shown [START_REF] Ladd | Lattice-Boltzmann Simulations of Particle-Fluid Suspensions[END_REF]] that this artefact does not induce consistent imprecisions in the simulation.

The hydrodynamic force applied on microspheres can also be computed based on the momentum exchange. It is indeed defined as:

F p = i 2f i (x, t + ) -6w i ρv p • e i • e i , (6.3)
where i denotes every boundary link belonging to the particle p.

Other slightly different definitions of equations (6.2) and (6.3) exist, for example in [START_REF] Feng | Coupled lattice Boltzmann method and discrete element modelling of particle transport in turbulent fluid flows : Computational issues[END_REF], [START_REF] Owen | An efficient framework for fluid -structure interaction using the lattice Boltzmann method and immersed moving boundaries[END_REF] and [START_REF] Wang | Evaluation of three lattice Boltzmann models for particulate flows[END_REF]. The differences concern from [START_REF] Owen | An efficient framework for fluid -structure interaction using the lattice Boltzmann method and immersed moving boundaries[END_REF] in particular the choice of the parameter that we define here equal to 6. This choice of parameter is taken from [START_REF] Leonardi | Coupled DEM-LBM method for the free-surface simulation of heterogeneous suspensions[END_REF] and generates the more accurate results in our computation of particle velocity.

Interaction between microspheres and between microspheres and wall As outlined previously in this section and in section 6.1, the scale of the phenomenon described here suggests that inter-particles collision forces do not play a fundamental role. Nevertheless, in order to present a more precise simulation, we decided not to neglect them.

As suggested in [START_REF] Feng | The immersed boundary-lattice Boltzmann method for solving fluid-particles interaction problems[END_REF], [START_REF] Feng | Coupled lattice Boltzmann method and discrete element modelling of particle transport in turbulent fluid flows : Computational issues[END_REF], such interaction can be considered in the same way as the interaction between microspheres and vessel wall.

Many contact laws are defined in the literature, some of the principal ones being the Hertzian contact law [START_REF] Leonardi | Coupled DEM-LBM method for the free-surface simulation of heterogeneous suspensions[END_REF][START_REF] Han | A combined finite/discrete element simulation of shot peening processes -Part I: studies on 2D interaction laws[END_REF] and the linear contact law [START_REF] Owen | An efficient framework for fluid -structure interaction using the lattice Boltzmann method and immersed moving boundaries[END_REF][START_REF] Feng | Coupled lattice Boltzmann method and discrete element modelling of particle transport in turbulent fluid flows : Computational issues[END_REF].

In this work, we decided to implement the method used by [START_REF] Basciano | Computer modeling of controlled microsphere release and targeting in a representative hepatic artery system[END_REF], which simulates such interaction through the definition of the restitution coefficient C R , defined by v + p = -C R v - p , where v + p and v - p are, respectively, the velocity of the particle p after and before the collision. [START_REF] Basciano | Computer modeling of controlled microsphere release and targeting in a representative hepatic artery system[END_REF] propose an analysis of microspheres transport 184CHAPTER 6. PATIENT-SPECIFIC MODELLING OF MICROSPHERES TRANSPORT simulation, testing different restitution coefficients for wall collision, showing that an ideal elastic collision, i.e. the choice of C R = 1, is reasonable.

Given two microspheres m and n (or a microsphere and the wall), we define from [START_REF] Hyakutake | Lattice Boltzmann simulation of blood cell behavior at microvascular bifurcations[END_REF] the impulsive force (the integral on the time step ∆t = 1 of the temporary collision force) of the colliding particle m as: (6.4) where M represents the particle mass, v rel the relative velocity between the two spheres (the velocity of the microsphere in the case of microsphere-wall contact), and n c is the normal unit vector, going from m to n.

F c = 1 2 (1 + C R )M v rel • n c n c ,
Newton's second law. Once the hydrodynamic force applied by the fluid on the particle is computed, the update of its position during time is governed by the Newton's second law:

ma + cv = F p + F c . (6.5)
The aim of the term cv is to take into account all the potential dissipation forces, representing a viscous force. In equation ( 6.5) we estimated the gravitational force to be negligible in this system.

In order to compute the particle velocity we integrate equation (6.5) using a central difference scheme, which is shown to be more stable than classical explicit Euler method.

It is shown in the literature by [START_REF] Feng | Coupled lattice Boltzmann method and discrete element modelling of particle transport in turbulent fluid flows : Computational issues[END_REF], [START_REF] Owen | An efficient framework for fluid -structure interaction using the lattice Boltzmann method and immersed moving boundaries[END_REF] among others, that a subcycling is needed for the integration of Newton's second law in the LBM scheme.

Indeed, two different time steps exist for the LBM scheme (∆t) and for the Newon's second law integration(∆t N ). In order to obtain a coherent simulation, the number N s of integration steps of equation (6.5), should verify ∆t = N s ∆t N . The behaviour of the numerical simulation as a function of N s is not clearly understood, but according to [START_REF] Feng | Coupled lattice Boltzmann method and discrete element modelling of particle transport in turbulent fluid flows : Computational issues[END_REF] and to our experience, the simulation is more stable when small values of N s are chosen.

Instabilities may occur when the particles velocity is high. Indeed, hydrodynamic and collision forces are not updated during the integration subcycling, and the particle is remapped on the lattice nodes only at the end of every LBM cycle. This suggests that, during one subcycle, microsphere position should not cross more than one lattice step.

It should be noted that we decided not to consider the lubrication force defined for instance in [Ladd, 1994a,b, Leonardi et al., 2014]. The choice is justified by the dilute nature

LATTICE-BOLTZMANN METHOD FOR THE SIMULATION OF PARTICLES TRANSPORT185

of the microspheres suspension in the saline solution. Such characteristic also allowed us to neglect the lubrication forces which can have an important role if the gap between two particles become very small, in order to correctly simulate the fluid compressed in between. Indeed, the lubrication force must be considered in the case of highly concentrated suspensions, i.e. if the solid fraction φ is greater than 0.5 [Ladd andVerberg, 2001, Stickel and[START_REF] Stickel | Fluid Mechanics and Rheology of Dense Suspensions[END_REF], which, as shown above in this section, is not the case.

Results

The method described above was tested on artificially defined geometries. In order to consider a condition similar to microspheres transport in small vessels, which requires a two-way simulation (i.e. considering the influence of microspheres on blood flow), spheres with a larger diameter than reality were considered. In particular, spatial discretisation h was defined as 0. ratio, between the first vessel and its descendants, appears in the 2D and 3D simulations visualisations. This is not surprising: in the 2D case, the velocity difference between the entry and the outlets is proportional to the ratio of the radii; in the 3D case it is proportional to the ratio between the squares of the radii. It was verified that, after a temporary insta- Microspheres tracking is achieved by solving the system (6.6) with a trapezoidal scheme.

The drag force is the only force considered here, neglecting gravity due to the very small particles size. The density of the glass microspheres was defined as 3600 Kg/m 3 , and their average diameter as 25 µm.

Results

Particles were injected at the mesh nodes of the hepatic artery input surface. The injection velocity was defined as 30 cm/s, equal to the input blood velocity. Indeed, Theraspheres ®

vendor suggests an injection modality assuring a velocity equal or superior to the blood velocity at the location where the catheter tip is placed. The trajectories of the particles injected in the hepatic artery for the four considered patients are shown in Fig. 6.6. It

Discussion and conclusion

Two different approaches for the simulation of microspheres trajectories are proposed.

Concerning the LBM, a two way model is proposed, taking into account the influence of the fluid on the microspheres trajectories and vice versa. The collision between particles is also considered. Due to its easy handling of complex structure and its potentially very efficient parallelisation, the contribution of this kind of model appears particularly interesting, under the condition that the particles size is non negligible relatively to the vessel size.

In such vessels, the LBM simulation appears very promising. The bases for such kind of future simulations are placed in this chapter.

On the other side, a numerically efficient one-way simulation of microspheres trajectories has been proposed using the FVM in the principal, segmented vessels. Indeed, in this case the influence of particles on blood flow can be neglected. This provides the amount of microspheres reaching every outlet of such principal tree. Since the principal vessels irrigating the tumours are identified, an estimation of the amount of microspheres reaching the tumour can be done. Nevertheless, this dramatically depends on the definition of the boundary conditions at the outputs of such tree. This is why the coupling of the presented simulation to the in vivo velocity measures acquired with MRI, and to the envisaged pressure measures, will be a priority for us. This will be complementary to the modelling of the hydraulic loads at the outlets of the principal tree, computed through the vascular model presented in Chapter 4 or using a lumped parameter model. In this context, the influence of the injection position on the spheres distribution will be analysed by comparing different injection sites, injection velocities and modalities. We will also analyse the influence, on the blood streamlines, of the catheter shape, distal position and angle relatively to the vessel direction, which is difficult to control exactly.

A preliminary validation of the LBM particles transport model was performed through a comparison with Ansys Fluent. The injection of few spheres of a diameter d s = 1.3 mm in a single, straight tube, of a dimension similar to the one of Fig. 6.3, was performed using the LBM and FVM proposed in the sections above. It was verified that, in both cases, the particles velocities reached a constant value corresponding to the velocity of the surrounding fluid.

The validation of the microspheres distribution represents an important challenge. During the injection of Theraspheres ® , no contrast product is added and their trajectories can not be visualised in any way. Their distribution is estimated based on the bremsstrahlung SPECT, whose accuracy is limited, as already mentioned in Section 2.1.3. Therefore, a possible validation could only be performed on the final distribution of microspheres in the entire tissue. Thus, we envisage to couple the present simulation to the entire simulated tree presented in Chapter 4, in order to estimate a dose distribution map in liver tissues.

Another perspective of the research team consists in the simulation of microspheres transport using a physical model. A 3D printing technique similar to stereolithography is envisaged, starting from the numerical tree model provided by the vascular growth model, or by the more realistic, but more limited, geometry introduced in Ansys Fluent.

Conclusion and perspectives

In this manuscript, the way to a computational patient-specific modelling of the SIRT treatment is opened. In order to dedicate particular attention to the clinical groundwork of the model, most of the different steps required for the development of a full model have been deeply analysed or formalised and tested under simplified conditions.

The clinical protocol was analysed, and improvements were proposed for an optimisation of the extraction of patient's data, capital for the numerical model. MRI sequences were tailored to the extraction of flow quantified data in the investigated vessels. The hepatic arterial tree is a particularly complex structure, due to its small size, rapidly decreasing along the first bifurcations, and due to its breathing-related important movement.

Even if the reliability of the proposed measurement technique needs to be confirmed on more patients, it still represents an important basis for the patient-specific computation of blood flow in the hepatic artery. In the future, the velocity data will be coupled to in vivo pressure measures, that we plan to acquire thanks to a specific catheter, placed along the path to the injection site.

The extraction of patient-specific quantified flow data, namely pressure and velocity, represents a crucial starting point for a numerical simulation of microspheres that attempts to be patient-specific. Indeed, the principal issue addressed in the literature, concerning numerical simulations in physiological domain, is the definition of the boundary conditions, associated to the partial derivative equations describing the evolution of the fluid. Such conditions have a direct impact on the simulation results, and at the same time they are often difficult or impossible to measure in vivo. The small size of the vessels considered in this thesis makes the task even more difficult, compared to most of the analogous studies, where bigger structures, typically in the cardiovascular context, are modelled.

As shown in Fig. 6.7, different modelling approaches were considered in this work, leading to preliminary or more advanced results, at each step of both approaches. On the one hand (left side of Fig. 6.7), a partial problem has been investigated, where the computational domain was considered as the principal hepatic arterial tree, down to a diameter permitting a reliable segmentation. Using the Finite Volumes Method and the Lattice Boltzmann Method, a blood flow simulation was achieved in such principal geometry, extracted from the patient's CT acquisitions. In the present work, the boundary conditions were not defined based on the patients data, but on values taken from the literature. The fluid simulation results obtained by the two methods were compared, leading to encouraging results that were coherent with the literature data.

Concerning microspheres transport in the patient's arteries, this was simulated using the FVM through the standard commercial software Ansys Fluent. The presented work will need to be completed, in order to analyse the influence, on the microspheres distribution, of different parameters as vessels geometry, blood pressure and velocity, as well as injection modalities. Considering the LBM, a method for a two way simulation of microspheres, considering inter-spheres collision, was developed. Encouraging preliminary results were obtained on a standard geometry, and with a simplified configuration of microspheres (in number and size).

For future patients, the proposed technique for the extraction of flow velocity will be used for the initialisation and the validation of the numerical simulation, leading to a full patient-specific simulation. The proposed procedure for the simulation of blood flow in the segmented arteries is currently limited by segmentation inaccuracies, which are not surprising due to the difficulties in accessing such vessel. In the future, beside an optimisation of the segmentation procedure, we plan to perform a sensitivity study in order to precise the influence of the vessels geometry, and in particular of the segmentation accuracy, on the flow profile. The injection step will also be analysed more deeply, in order to investigate the influence of the catheter presence on the fluid streamlines and the possible influence of the saline solution injected with the microspheres on blood viscosity. This will let us estimate the error inherited in the simulation and strengthen the reliability of our simulations, quantifying inaccuracy of microspheres distribution.

As far as validation is concerned, the two CFD methods will have to be applied on more patients and the results will have to be compared. Moreover, simulation results will be compared to measures acquired on a resin printed physical model.

On the other hand, a more complete problem has been investigated (cf. right side of Fig. 6.7). The vascular growth model was accommodated to a patient-specific simulation of the patient's complete vascular tree, based on its anatomical characteristics. The importance of the definition of BC is central again: input and output pressure and flow values have a direct impact on the structural and functional characteristics of the complete, simulated tree. Thanks to multi-GPU parallelisation, such tree can be generated down to the arterioles level and might then be used as a computational domain for a LBM simulation of blood flow and microspheres transport in vessels much smaller than the segmented ones. Such simulations will then be validated by comparison of the obtained microspheres distribution until the smallest vessels, with the real distribution given by post-treatment scintigraphy.

Since approximate blood pressure and flow values are available at every bifurcation of the simulated complete tree, this could provide more realistic boundary conditions at the outputs of the principal tree for 3D CFD simulation.

Concerning CFD modelling of the complete arterial tree, another perspective is to The relationship between the strength of the main magnetic field and the frequency of precession ω 0 is given by the following equation, called the Larmor equation:

ω 0 = γ → B 0 , (A.1)
where γ corresponds to the gyromagnetic ratio associated to the hydrogen proton and ω 0 is called the Larmor frequency.

Radiofrequency energy

In order to acquire signal from the tissues, the equilibrium state described in the previous paragraph is excited by a radiofrequency (RF) wave sent towards the patient (Fig. A.1, right). The mechanism is the following. The energy transmission from the RF pulse to the spinning protons is maximal if the radiofrequency wave is at the same precessional frequency than hydrogen protons, the Larmor frequency. In this manner, resonance occurs between the RF pulse and the spinning protons that start precessing all in phase, creating a non-zero transverse magnetisation → M xy . When the RF energy, sent perpendicularly to

T1 and T2 relaxation times and contrasts

The signal is acquired while the magnetisation returns back to the equilibrium phase.

This phenomenon is called "relaxation" and it is the moment when every tissue behaves in different manners and can be distinguished. Moreover, the moment where the signal is acquired during the recovery of the equilibrium determines the type of contrast of the acquired image. The signal detected from the MR system is called Free Induction Decay The particularity of the momentum equation is that the transported scalar quantity is a velocity component, and the velocity vector is unknown itself. Moreover, the source term depends on the pressure field, which is also unknown. In order to solve this hindrance, pressure-velocity coupling algorithms are used.

First, continuity and momentum equation are written in a discretised form. The same discretisation scheme leading to (B.14) can be applied to the momentum equations. A central differences scheme is used for the discretisation of the diffusion term in order to write it in terms of the values of the velocity vector at the centres of the surrounding cells, u L , L ∈ N (K). The advection term can be discretised using several schemes, common ones being the first or second order upwind scheme. The discretised equation for the x momentum equation can thus be written as: (B.15) where the subscript f refers to the cell face f , A f is the d-dimensional generalisation of the area of the face f and N f is the number of faces of the cell K and corresponds, except on the boundaries of the domain, to the number of elements in N (K) . The quantity p f represents the value of p on the cell face f and its expression in terms of the surrounding values of pressure can be achieved according to different methods. The simplest is a linear method, where the p f is computed as the average of the values of pressure at the centres of the neighbouring cells. The default scheme in Ansys Fluent is the second order scheme, based on a central differences scheme. Face pressure is computed as: (B.16) where p L is the pressure evaluated at the centre of the cell L. for the face K|L is computed as:

a p u p = L∈N (K) a L u L + N f f p f A f • î + S,
p f = 1 N f L∈N (K) p L + 1 N f L∈N (K) ∇p L • -→ r L ,
J f = ρ a p,K v n,K + a p,L v n,L a p,K + a p,L + d f (p K + ∇p K • -→ r K ) -(p L + ∇p L • -→ r L ) (B.18) = Ĵf + d f p K -p L ,
where v n,K , v n,L are the normal velocities within cells K and L and d f a coefficient taking into account the average of the coefficients a p of the momentum equations for the cells K and L. This expression is finally replaced in the continuity equation.

In order to solve the system of the Navier-Stokes equations, an iterative algorithm is used, a very common one being the SIMPLE algorithm original developed by [START_REF] Patankar | A Calculation Procedure for Heat, Mass and Momentum Transfer in Three Dimensional Parabolic Flows[END_REF]. The SIMPLE algorithm is a predictor-corrector method, where the initial guess of pressure and velocity are corrected based on the continuity equation. Thanks to the procedure described above, the continuity equation includes pressure and velocity, and can thus be used to correct the pressure field. Indeed, the velocity computed from the momentum equation using the correct pressure field should satisfy the continuity equation.

Otherwise, the continuity equation is used as a pressure correction equation. The different steps of the algorithm are presented below:

1. A guessed pressure field p * is used to compute, through the momentum equation, a temporary velocity field u * .

2. Using the relationship (B.18), the temporary flux J * f is computed, but at this step it does not verify continuity. It is thus corrected by J f = J * f + J f , which satisfies continuity. Defining as p the cell pressure correction, then the flux correction is written in the SIMPLE scheme in the form:

J f = d f p K -p L .
3. The newly defined flux correction equation is substituted in the discretised continuity equation (B.17) in order to obtain the pressure correction equation for p :

a p p = L∈N (K) a K p K + b, (B.19)
where the source b corresponds to the net flow rate into the cell: b = N f f J * f A f . Such pressure correction equation is solved using the algebraic multigrid method [START_REF] Hutchinson | A multigrid method based on the additive correction strategy[END_REF].

J f = J * f + d f p K -p L , (B.21)
where α p is the under-relaxation factor, which reduces the change of pressure during each iteration.

5. If convergence is not reached, the previous steps are repeated.

Abstract

Selective internal radiation therapy using Yttrium-90 loaded glass microspheres injected in the hepatic artery is an emerging, minimally invasive therapy of hepatocellular carcinoma, which is the fourth cause of mortality in the world. Currently, microspheres distribution can be only approximately predicted by the injection of a radiotracer, whose behaviour may be different. A personalised intervention can lead to high concentration dose in the tumour, while sparing the surrounding parenchyma. This work is concerned with the development of a patient-specific numerical model for the simulation of microspheres trajectories and treatment optimisation. Clinical imaging protocol is utilised and optimised in order to extract patient's specific data such as liver, tumours, hepatic artery and blood flow. Normal and malignant hepatic arterial vasculature and tissues are simulated down to a vessels diameter of 0.05 mm. A preliminary simulation of microspheres distribution in liver tissue is proposed and validated against post-treatment scintigraphy. Microspheres are here supposed to distribute proportionally to blood flow, which is computed based on Poiseuille's law. More precise computational fluid dynamics (CFD) simulations of blood flow in the patient's segmented arteries are performed. The Finite Volume Method (Ansys Fluent) and the Lattice Boltzmann Method (in-house developed software) are used to this purpose and their efficacy is compared. Microspheres transport is simulated in the patient's hepatic artery using the FVM, and in a representative geometry using the LBM method.

A phase contrast MRI sequence has been optimised in order to extract blood velocity from the hepatic artery and validate CFD simulations.

Résumé

La radiothérapie interne sélective est une thérapie émergente très peu invasive du carcinome hepatocellulaire, quatrième cause de décès par cancer dans le monde. Des millions de microsphères chargées en Yttrium 90 sont injectées dans l'artère hépatique par un cathéter.

  Merci également à Jean-Claude Nunes, qui m'a appris les bases du recalage d'images et offert son aide dans des moments critiques. Je voudrais ensuite remercier les rapporteurs Gabriele Dubini et Frédérique Frouin, ainsi que la présidente de mon jury de thèse Dominique Barthès-Biesel, pour le temps qu'ils ont dédié à une analyse approfondie et curieuse de ce manuscrit, pour leur commentaires constructifs et leur participation active et fondamentale lors de ma soutenance. Merci à Simon Esneault et à Gemil Goksu, président de l'entreprise Therenva, pour nous avoir permis d'utiliser leur logiciel Endosize et nous avoir aidés à l'utiliser lors du démarrage de ma thèse. Merci à Michel Rochette et à Valery Morgenthaler pour m'avoir permis d'utiliser le logiciel Ansys et pour tout leur soutien. Merci à Antoine, Maëva et Paul, les stagiaires qui m'ont accompagné pendant des périodes plus ou moins courtes. Merci pour vos travaux et merci à ceux qui ont eu la chance de faire cobaye dans l'IRM !

  multidisciplinary approach proposed in this dissertation for the optimisation of liver targeting in SIRT protocol . . . . . . . . . . . . . . 1.1 Views of the front, visceral and lower surfaces of the liver. [Virtual Liver: 3D Liver Anatomy, 2012] . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Liver is divided into eight segments by the hepatic veins [Hitachi Medical Systems America Inc., 2012] . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Liver tissue as an assembly of small hexagonal lobules . . . . . . . . . . . . 1.4 Major arteries sprouting from the celiac trunk [Elsevier Images, 2016] . . . 1.5 Diagram of velocity and pressure in human systemic circuit [Martini et al., 2012] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6 Injection of radioactive microspheres into the hepatic artery and their deposition in the tissues . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.7 After the injection of 99m Tc-MAA, a SPECT-CT is realised. Injected liver volume is measured as shown on the top image based on the activity distribution. Tumour volume is measured using a higher threshold (bottom image) [Garin et al., 2012] . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.8 Outline of the multidisciplinary steps needed for a complete simulation of SIRT treatment through an integrative approach with their main bibliographic references . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.9 Overview of the mathematical modelling and simulation of blood flow and microspheres transport considered in Chapters 4, 5 and 6 of this thesis . . . 2.1 Scheme comparing a 1D detector fan-beam X-rays CT-scan and 2D detector cone-beam X-rays CT-scan [Scarfe et al., 2006] . . . . . . . . . . . . . . . . LIST OF FIGURES 2.2 Three CT images acquired during the distribution of contrast agent in liver tissues and vasculature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Scheme of the detection of coincident photons from a PET gamma camera . 2.4 Comparison between a SPECT acquired after the injection of 99m Tc-MAA and a bremsstrahlung SPECT acquired after the 90 Y radiotherapy . . . . . 2.5 Example of the axial, sagittal and coronal views of bremsstrahlung SPECT acquired after the injection of 90 Y microspheres, CT and their superimposition. Images acquired at CEM . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 Two examples of angiographies acquired at CEM in the Theraspheres ® framework, guiding the choice of the optimal catheter placement for the injection of the 99m Tc-MAA . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 Example of dynamic gadolinium enhanced MRI realised in the case of a HCC 2.8 Scheme of the imaging protocol applied in the case of patients candidates for a Theraspheres ® treatment . . . . . . . . . . . . . . . . . . . . . . . . . 2.9 Examples of artefacts related to phase wrapping of two phase contrast magnitude image . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.10 The effect of the bipolar gradient on stationary spins and on moving spins during a phase contrast acquisition . . . . . . . . . . . . . . . . . . . . . . . 2.11 Diagram of a typical phase contrast sequence . . . . . . . . . . . . . . . . . 2.12 Aliasing on the vena cava . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.13 Repetition of "phantoms" of aorta and vena cava on a 3D PC sequence acquired with TE=4.08 ms, TR=18 ms and a slice thickness of 4 mm . . . . 2.14 Screenshots of the DICOM header of the 3D phase contrast sequence realised on the 1.5 T Optima MR450w GE system of the CEM centre . . . . . . . . 2.15 Reducing the receiver bandwidth leads to a better SNR, while increasing acquisition time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.16 3D PC sequence was tested on a plastic tube filled with flowing saline solution 2.17 MR images acquired with the optimised 3D PC sequence . . . . . . . . . . . 2.18 MR images acquired with the optimised 2D PC sequence . . . . . . . . . . . 2.19 Case of a patient with a standard geometry of hepatic arteries . . . . . . . . 2.20 For this patient, only the beginning of the common hepatic artery, sprouting from the celiac trunk, is visible . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Schematic representation of an arteriovenous shunt induced by the presence of an HCC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Example of ROI selection on a planar scintigraphy . . . . . . . . . . . . . . 3.3 Example of a CT-scan at arterial phase with an early enhancement of the hepatic vein . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Schematic representation of a coronal view of liver and the principal related veins . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Massive tumour invading vena cava . . . . . . . . . . . . . . . . . . . . . . . 3.6 Example of the correct placement of the three principal ROIs: on the hepatic vein, and on the vena cava, before and after the connection to hepatic veins 3.7 Example of the sequence of T1-weighted MRI acquisitions during the diffusion of contrast agent . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 Placement of the ROIs on the three hepatic veins (right, middle and left) on the three LAVA images acquired during the arterial phase. . . . . . . . . 3.9 The mean and standard deviation of the hepatic veins enhancement at the arterial time are plotted as a function of the shunt fraction . . . . . . . . . . 3.10 The ratio between the hepatic veins enhancements at the portal time is plotted based on the relative shunt fraction . . . . . . . . . . . . . . . . . . 3.11 The difference in signal intensity between the two ROIs placed on the vena cava before and after the arrival of hepatic veins is plotted based on the relative lung shunt fraction . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.12 Enhancement curves of the three hepatic veins (left LHV, median MHV and right RHV) during the arterial phase, realised through the analysis of the multi-arterial MRI acquisitions. . . . . . . . . . . . . . . . . . . . . . . . . . 3.13 CT image acquired during the arterial phase. The presence of an extensive tumour prevents the localisation of the mean hepatic vein. Other hepatic veins are marked by yellow ROIs . . . . . . . . . . . . . . . . . . . . . . . . 3.14 Comparison between two LAVA sequences tested on a healthy volunteer . . 4.1 Outline of the method used for a patient-specific simulation of the hepatic arterial tree . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Two slices of a MR image acquired at CEM during the portal phase, showing the principal challenges in liver segmentation. . . . . . . . . . . . . . . . . . LIST OF FIGURES 4.3 Surface evolution during four cycles of the liver shape segmentation in MESA software . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Example of binary tree ordered in a "depth-first" manner and of its associated input file . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Example of semi-manual segmentation of vessels with EndoSize. At each step the user defines a new point inside the vessels. Unlike the visual representation proposed in this picture, the vessel centrelines are considered as piecewise straight lines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Outline of the workflow leading to the vascular model input file . . . . . . . 4.7 Example of different contrast enhancements in 3D angiography and CT-scan 4.8 Registration process between portal phase CT-scan and 3D Cone Beam angiography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.9 General overview of the model simulating the hepatic vascular trees and tissue131 4.10 Flow chart representing two loops of events, distinguished in the simulation of the organ (tissue and vascular network) growth . . . . . . . . . . . . . . . 4.11 Perfusion process of newly created macro-cell, by hepatic vascular network . 4.12 The binary vascular tree, made of successive bifurcations . . . . . . . . . . . 4.13 Method used to ensure the consistency of a bifurcation . . . . . . . . . . . . 4.14 3D Cone-Beam angiography (left) and arteries segmentation in EndoSize (right) for both patients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.15 Example of a binary tree automatically extracted from the arteries segmented with the Frangi's filter on Patient 1 . . . . . . . . . . . . . . . . . . 4.16 Results of liver segmentation in MESA: final 3D surface (upper left) and 2D result (lower right) for 2 patients . . . . . . . . . . . . . . . . . . . . . . . . 4.17 Registration of arterial CT-scan (right) on portal CT-scan (left). Checkerboard after registration is shown below . . . . . . . . . . . . . . . . . . . . . 4.18 Registration of angiography on arterial phase CT-scan is determined from Frangi's enhanced 3D volumes. Superimposition of the two volumes is shown below the two images . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.19 Registration results between CT-scan and 3D angiography for both patients. Superimposition of liver and tumour contours obtained from portal phase CT and vascular skeleton segmented from angiography . . . . . . . . . . . . 4.20 Entire simulated hepatic arterial trees for both patients, including tumour vascularisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.21 Results of microspheres distribution simulation, based on Poiseuille's law. . 4.22 Comparison between simulated and actual microspheres distribution in liver tissue . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.23 Comparison between the automatically extracted arterial tree of Patient 1 and the larger simulated vessels . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 The two most common 2D and 3D LBM models for the velocity discretisation. Images from [Ehrhardt, 2013] . . . . . . . . . . . . . . . . . . . . . . . 5.2 Examples of extraction of the principal vessels of the simulated tree, according to different diameter criteria . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 LBM blood simulation in the patient-specific arterial tree of Patient 1 . . . 5.4 LBM blood simulation in the patient-specific arterial tree of Patient 2 . . . 5.5 An example of the procedure developed for the identification of the input and output surfaces as independent meshes . . . . . . . . . . . . . . . . . . 5.6 Cone-beam angiography and blood flow simulation results on segmented arterial tree for Patient 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.7 Cone-beam angiography and blood flow simulation results on segmented arterial tree for Patient 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.8 Cone-beam angiography and blood flow simulation results on segmented arterial tree for Patient 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.9 Cone-beam angiography and blood flow simulation results on segmented arterial tree for Patient 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.10 Comparison between velocity fields obtained with the FVM and the LBM in the extracted hepatic arterial tree of Patient 1 . . . . . . . . . . . . . . . 5.11 Comparison between velocity fields obtained with the FVM and the LBM in the extracted hepatic arterial tree of Patient 2 . . . . . . . . . . . . . . . 5.12 Example of complex vascular structures where the parallelised LBM simulation could simulate blood flow. The diameter of the smallest vessels in the presented geometry is 0.5 mm . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1 Lattice representation of a 2D sphere . . . . . . . . . . . . . . . . . . . . . . List of Tables 2.1 Parameters of the 3D phase contrast sequence used to measure fluid velocity on the test object . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Parameters of the optimised 2D and 3D phase contrast sequences we propose for the quantification of blood velocity in the hepatic artery . . . . . . . . . 2.3 Blood velocity measured on the proposed optimised 3D PC MRI sequence on a healthy volunteer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Blood velocity measured on the proposed optimised 2D PC MRI sequence on a healthy volunteer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Velocity measurements realised on the 3D and 2D PC optimised sequenced acquired on 8 HCC patients (P1, ..., P8) at CEM and the average value. HA stands for hepatic artery, PV for portal vein and VC for vena cava. Values are given in cm/s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 Comparison of average blood velocities in aorta and vena cava measured with US, 2D PC MRI and 3D PC MRI on the same volunteer . . . . . . . .

  coefficient and the parameters of linear regression approximation concerning all the values measured on the CT-scan data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Results of the F-test and Student's t-test . . . . . . . . . . . . . . . . . . . . LIST OF TABLES 3.4 Confusion matrix obtained for the validation of the classification of the population according to the threshold given by a shunt percentage of 20%. Every patient is classified once: the sum of the values given in the matrix represents the number of patients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 3.5 Parameters of the rapid multi-arterial LAVA sequence for shunt detection, producing 4 clustered images, every 9 seconds, on the coronal plane . . . . . 117 4.1 Geometry of the CT acquisitions used for segmentation of hepatic arteries and liver shape. Values are given in mm . . . . . . . . . . . . . . . . . . . . 136 4.2 Frangi's filter parameters for hepatic arteries enhancement on cone-beam angiography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 4.3 Frangi's filter parameters for hepatic arteries enhancement on liver arterial CT-scan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 4.4 Input parameters for the simulation of the hepatic arterial tree. Initial number of macro-cells corresponds to the number of leaves of the segmented arterial tree . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 4.5 Percentage of the injected activity which deposits to the tumour for both patients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 5.1 Mesh characteristics for the hepatic arteries of the four patients . . . . . . . 168 5.2 Radius, area-weighted average pressure and velocity results in segmented arterial trees of the 4 patients . . . . . . . . . . . . . . . . . . . . . . . . . . 170 5.3 Blood flow velocity measured on phase contrast MRI acquisition (pixel size 0.9 mm) on a 27 years old healthy volunteer . . . . . . . . . . . . . . . . . . 173 Résumé de la thèse 1 Introduction et contexte médical

  La radiothérapie présente des contre-indications importantes à cause des effets secondaires des irradiations dans les tissus environnants. La chimiothérapie est limitée par le développement d'une résistance au traitement et par de nombreux effets secondaires. Des traitements locaux intra-artériels ont été développés dans le but d'améliorer le ciblage de la tumeur et réduire les effets secondaires. Ce type de traitements exploitent la structure particulière de la vascularisation hépatique: le foie est irrigué environ à 75% par la veine porte et à 25 % par l'artère hépatique. Le cancer du foie est irrigué principalement par l'artère hépatique et très peu par la veine porte. De plus, le phénomène d'angiogenèse tumorale fait que la vascularisation des régions cancéreuses est beaucoup plus importante que dans le parenchyme sain. Il est alors possible d'accéder à la zone tumorale en passant par le réseau artériel avec un cathéter introduit au niveau de l'artère fémorale. Grâce à l'apport important du sang à la tumeur, le traitement a tendance à être capturé plutôt par la tumeur que par le tissu sain. Des traitements d'embolisation simple des vaisseaux irriguant la tumeur, de chimioembolisation ou de radioembolisation sont ainsi possibles. L'administration d'un agent chimiothérapeutique directement dans les artères irriguant la tumeur est maintenant considérée comme la référence pour le traitement du cancer du foie à un stade intermédiaire, souvent couplé à une embolisation des vaisseaux irriguant la tumeur.

  agrégats d'albumine (MAA) chargés au Technétium 99m ( 99m Tc), un radiomarqueur ayant pour but d'imiter le comportement des Theraspheres ® et d'en prévoir la distribution. Une scintigraphie SPECT/CT (en anglais single photon emission computed tomography) est ensuite réalisée pour déterminer et quantifier la fixation du matériel radioactif. Cette étape permet de vérifier l'absence de la principale contre-indication au traitement, c'est à dire une fixation extra-hépatique, et de déterminer la captation de la tumeur afin de prévoir la dose d'Yttrium à injecter. En absence de contre-indication, il est possible de réaliser la radioembolisation à l' 90 Y. Le radiologue replace le cathéter à l'endroit choisi lors de l'artériographie diagnostique et le médecin nucléaire injecte les Theraspheres ® . Une autre scintigraphie est réalisée une heure après pour connaître l'étendue du foie effectivement traité, et des images scanner ou IRM sont réalisées dans les mois suivants dans le but de contrôler l'évolution de la tumeur. Le protocole utilisé actuellement présente de larges possibilités d'optimisation, concernant surtout la personnalisation du traitement, ainsi que le caractère invasif de la réalisation de deux angiographies. De plus, le choix du point d'injection est actuellement fait de façon essentiellement qualitative, et dépend de l'expérience du radiologue. Dans cette thèse une approche basée sur l'imagerie clinique a été adoptée dans le but de développer un modèle numérique de simulation du traitement par RTMI, afin d'adapter la distribution des microsphères à chaque patient. Une telle approche multidisciplinaire doit s'intéresser d'abord à l'extraction de données des patients à partir des images cliniques.

  et est totalement non invasive (sauf dans les rares cas d'allergie au produit de contraste). Cette dernière présente aussi de nombreux avantages grâce à sa meilleur flexibilité dans le choix des plans d'acquisition et à la possibilité d'acquérir plus d'images pendant la distribution du produit de contraste. L'utilisation de l'IRM est encore limitée par un coût important et par des durées d'acquisition parfois très longues. Un travail d'optimisation d'une séquence d'IRM de contraste de phase a été mené dans le but de l'adapter à la mesure de la vitesse du sang dans l'artère hépatique. Le flux sanguin dans l'artère hépatique est très peu connu et seules des informations approximatives de vitesse et de pression à l'entrée de celle-ci sont disponibles dans la littérature. L'acquisition d'information de vitesse le plus loin possible dans l'arbre artériel a un intérêt capital dans le développement de simulations spécifiques-patient de l'écoulement du sang. En effet, une simulation est considérée comme réaliste seulement si les bonnes conditions aux limites (en particulier en entrée et en sortie) ont été utilisées. Une vérité terrain concernant la vitesse et/ou la pression du sang à différents endroits permettrait une validation des conditions aux limites choisies.

Figure 1 :

 1 Figure 1: Outline of the multidisciplinary approach proposed in this dissertation for the optimisation of liver targeting in SIRT protocol

Figure 1 . 1 :

 11 Figure 1.1: Views of the front, visceral and lower surfaces of the liver. [Virtual Liver: 3D Liver Anatomy, 2012]

Figure 1 . 2 :

 12 Figure 1.2: Liver is divided into eight segments by the hepatic veins [Hitachi Medical SystemsAmerica Inc., 2012]

Figure 1 . 3 :

 13 Figure1.3: Liver tissue is an assembly of small hexagonal lobules whose size varies between 1 and 2.5 mm. The centre of every lobule is occupied by a central vein. Every corner of the hexagon corresponds to a portal triad, made of an arteriole, a portal venule and a bile duct[START_REF] Patton | Anatomy and physiology[END_REF]

Figure 1 . 4 :

 14 Figure 1.4: Major arteries sprouting from the celiac trunk [Elsevier Images, 2016]

Figure 1 . 5 :

 15 Figure 1.5: Diagram of velocity and pressure in human systemic circuit[START_REF] Martini | Fundamentals of Anatomy and Physiology[END_REF]

Figure 1 . 6 :

 16 Figure 1.6: Injection of radioactive microspheres into the hepatic artery (left, [BTG IM -Therasphere RoW, 2016]) and their deposition in the tissues (black circles), the bar representing 1 mm (right,[START_REF] Kennedy | Pathologic response and microdosimetry of (90)Y microspheres in man: review of four explanted whole livers[END_REF]

FIGURE 1 .

 1 FIGURE 1. Delineation of VOIs used for quantitative analysis of SPECT/CT data. VOI defined on the injected liver (A) and tumor (B).

FIGURE 1 .

 1 FIGURE 1. Delineation of VOIs used for quantitative analysis of SPECT/CT data. VOI defined on the injected liver (A) and tumor (B).

Figure 1 . 8 :

 18 Figure 1.8: Outline of the multidisciplinary steps needed for a complete simulation of SIRT treatment through an integrative approach with their main bibliographic references

Figure 1 . 9 :

 19 Figure 1.9: Overview of the mathematical modelling and simulation of blood flow and microspheres transport considered in Chapters 4, 5 and 6 of this thesis.

Figure 2 . 1 :

 21 Figure 2.1: Scheme comparing a 1D detector fan-beam X-rays CT-scan and 2D detector conebeam X-rays CT-scan[START_REF] Scarfe | Clinical Applications of Cone-Beam Computed Tomography in Dental Practice[END_REF]

Figure 2 . 2 :

 22 Figure 2.2: CT images acquired during the distribution of contrast agent in the liver. (a) Before the injection, aorta (1) and the tumours (2, 3) have lower signal intensities than hepatic parenchyma. (b) During the arterial phase aorta and hepatic arteries are enhanced, as well as the tumours vascularised by the hepatic artery. (c) During the portal phase contrast agent reaches the liver through the portal vein, which is enhanced as well as all the hepatic parenchyma. An extended necrotic region can be noticed inside the larger tumour (2). Images acquired at CEM, Rennes

Figure 2 . 3 :

 23 Figure 2.3: The coincident pair of photons is detected by the collimator of the gamma camera along the line or response (LOR). The photons are considered coincident if they are detected within a timing window of length τ .

  photons moving away in opposite directions along the same line, called the line of response (Fig. 2.3). The trajectory of the coincident pair of photons is detected by the collimator of the gamma camera, which is able to deduce the 3D position of the annihilation and thus the concentration of the radionuclide in every point of the tissue.The physical principles of SPECT are very similar, except for the used radiotracer. In this case the radioactive isotopes emit directly gamma rays, which are detected by the gamma camera that rotates around the

  Fig. 2.5 gives an example of the axial, sagittal and coronal views of 90 Y bremsstrahlung SPECT, CT and their superimposition. SPECT and PET are commonly called scintigraphy, the name being related to the employed scintillation camera.

  Figure 2.6: Two examples of angiographies acquired at CEM in the Theraspheres ® framework, guiding the choice of the optimal catheter placement for the injection of the 99m Tc-MAA

Figure 2 . 7 :

 27 Figure 2.7: Example of dynamic gadolinium enhanced MRI realised in the case of a HCC. Soon after the injection, arteries and arterially-perfused tissues like HCC are enhanced. During the portal phase the maximal enhancement of hepatic parenchyma is observed. Later, contrast agent filters into the extracellular space

Figure 2 . 8 :

 28 Figure 2.8: Scheme of the imaging protocol applied in the case of patients candidates for a Theraspheres ® treatment

Figure 2 . 9 :

 29 Figure 2.9: Examples of artefacts related to phase wrapping of a sagittal (left image) and a axial (right) phase contrast magnitude images. The field of view does not cover the whole body in the phase encoding direction (anterior to posterior in the left image, right to left in right one), thus the excluded posterior (left image) and lateral (right image) parts are superimposed on the image.In the right image, arms can clearly be recognised

Figure 2 .

 2 Figure 2.10: The effect of the bipolar gradient on stationary spins and on moving spins during a phase contrast acquisition. The application of the second lobe of the bipolar gradient (right image) cancels out the phase shift of the stationary spins while it increases the one of moving protons.

Figure 2 .

 2 Figure 2.11: Diagram of a typical phase contrast sequence. The two opposite lobes of a bipolar gradient (dotted lines) can be applied on any of the three directions. Simultaneous application of different gradients can lead to encode velocity on any direction.

Figure 2 .

 2 Figure 2.12: Aliasing on the vena cava. The diagram shows the grey value curve calculated on the red line. Venous flow has a laminar character, and velocity in the centre is higher. The central pixels are darker because here the velocity exceeds Venc

Figure 2 .

 2 Figure 2.13: Repetition of "phantoms" of aorta and vena cava on a 3D PC sequence acquired with TE=4.08 ms, TR=18 ms and a slice thickness of 4 mm

Figure 2 . 14 :

 214 Figure 2.14: Screenshots of the DICOM header of the 3D phase contrast sequence realised on the 1.5 T Optima MR450w GE system of the CEM centre

Figure 2 . 15 :

 215 Figure 2.15: Reducing the receiver bandwidth ∆F leads to a better SNR (image noise grows as √ ∆F ), while increasing acquisition time. If the FOV is to be unchanged, the number of samples is constant and a smaller bandwidth induces a longer sampling time ∆t, according to the relationship ∆F = 1 ∆t

( a )

 a Figure 2.16: 3D PC sequences were first tested on a plastic tube filled with flowing saline solution, in order to validate velocity measurements. On the right image, The greyscale gives velocity values on the x direction. Dark pixels show velocity going from the reader to the text, white pixels mark velocity going towards the reader

  2.2. PROPOSED IMPROVEMENTS OF THE CURRENT CLINICAL PROTOCOL 87 and Fig. 2.18. The reconstructed velocity maps (cf. equation (2.1)) are shown as well. Concerning the 3D sequence, we decided to output only the magnitude image and the one offering the absolute velocity in every direction in order to reduce the post-processing time of the 3D sequence during the acquisition of the clinical protocol. In the 3D image (a) Magnitude-weighted phase image (b) Magnitude image (c) Reconstructed velocity map

Figure 2 .

 2 Figure2.17: MR images acquired on healthy volunteer using the proposed optimised 3D PC sequence whose parameters are given inTable 2.2

  Figure 2.18: MR images acquired on healthy volunteer using the proposed optimised 2D PC sagittal sequence whose parameters are given in Table 2.2. In (a), (b) and (c) the contrast is proportional to the velocity in the given direction. White pixels indicate velocity in the positive direction, black pixels in the negative one. In (d) is the phase image proposed by the MR system and (e) is the magnitude image. (f ) represents the reconstructed image giving a pixel value equivalent to velocity in mm/s. The yellow arrow points to the hepatic artery on all the six images, the green one to the test object

  Figure 2.19: Case of a patient with a standard geometry of hepatic arteries. (a) Common (CHA), proper (PHA), left (LHA) and right (RHA) hepatic arteries are well visible on the same slice. The sprouting of the gastroduodenal artery (GdA), which separates the common hepatic artery to the proper one is also visible. (b) MIP image of all the acquired region: vena cava (VC), aorta, portal vein (PV), one renal artery (RA) and an hepatic vein (HV) of the patient are clearly visible. A part of the common hepatic artery (CHA) can be seen next to the portal vein.

  Figure 2.20: For this patient (P2 in Table 2.5), only the beginning of the common hepatic artery (CHA), sprouting from the celiac trunk, is visible. Even on the MIP image (right), the visibility of the common hepatic artery decreases when leaving the aorta. Vena cava (VC), aorta, renal arteries (RA) and splenic artery (SA) of the patient are clearly visible

  A posteriori analysis of the clinical images of patients diagnosed with SIRT contraindication factors can indeed help in the establishment of indicative parameters and of new imaging methods for the earlier diagnosis of such factors.Side effects of SIRT are numerous, and can be classified into three categories: complications related to the surgical procedure, radiations to healthy hepatic parenchyma and extra hepatic microspheres deposition.A surgical operation always presents certain risks, for example allergic reactions to contrast product. Nevertheless, vascular or infectious complications are very rare in a specialised centre. An important issue to be taken into account in SIRT concerns the vascular modifications that can be induced by the presence and the placement of catheter.Vessel stimulation can induce the development of spasms, blood clots that will flow into the vascular network, and modifications in the geometry of the vessels induced by the deployment of a rigid catheter guide. The preservation of the vessels geometry is clearly essential for the proper performance of the treatment, which is tailored based on the distri-CHAPTER 3. AT-RISK PATIENTS IDENTIFICATION bution of 99m Tc-MAA in the vascular network geometry at the moment of the diagnostic angiography.

  Fig. 1.4), and in the gallbladder. Radiations reaching the gastrointestinal tract can lead to the development of ulcers. If the diagnostic angiography reveals the presence of gastrointestinal deposition, radiologist has two choices. Realise a second angiography with a different catheter placement, or physically occlude the vessel involved. The latter option is clearlythe more invasive and is realised using surgical clips or through the deployment of a metallic coil in the vessel. The vessels that physicians decide to occlude are never terminal ones, meaning that other vessels can provide the needed blood amount to the organ. Radiations induced gallbladder inflammation can appear if microspheres deposition occurs through the cystic artery, which is a branch of the proper hepatic artery (cf. Fig.1.4). The cystic artery cannot be permanently embolised due to its terminal character: no other vessel contribute to gallbladder oxygenated blood supply. An alternative considers the stimulation with a catheter of the concerned vessel, in order to generate spasms that will temporary close the vessel.

Figure 3 . 1 :

 31 Figure 3.1: Schematic representation of an arteriovenous shunt induced by the presence of an HCC. The tumour (black ellipse) can create a link between the arterial and venous network due to its angiogenic character

Figure 3 . 2 :

 32 Figure 3.2: Example of ROI selection on the anterior and posterior projection of a planar scintigraphy acquired one hour after the injection of 99m Tc-MAA. Green and blue ROI define, respectively, right and left lungs, while the red one indicates the liver. The ROI defining the mediastinum is not represented here. The image was acquired at the Eugène Marquis centre

Figure 3 . 3 :

 33 Figure 3.3: Example of an arterial CT with early enhancement of the hepatic vein (ROI). The enhanced region of the liver around the hepatic vein marks the presence of a massive tumour. (1) Inferior vena cava is not enhanced at arterial phase. (2) High contrast uptake in aorta. (3) Vertebral muscle. (4) Vertebrae. (5) Spleen

Figure 3 . 4 :

 34 Figure 3.4: Schematic representation of a coronal view of liver and the principal related veins.Yellow ellipses mark the placement of the two ROIs on the vena cava, before and after the connection with hepatic veins, and the ROI placed on the hepatic vein.

Figure 3 . 5 :

 35 Figure 3.5: The picture represents an arterial phase CT-scan of a massive tumour (marked in red on the right image) invading the vena cava (blue on the right image). In such circumstances, it is almost impossible to identify the vena cava when it crosses the liver

Figure 3 . 6 :

 36 Figure 3.6: Example of the correct placement of the three principal ROIs: on the hepatic vein, and on the vena cava, before and after the connection to hepatic veins

Figure 3 . 7 :

 37 Figure 3.7: Example of the sequence of T1-weighted MRI acquisitions during the diffusion of contrast agent. A LAVA-Flex sequence is performed before the injection (a). The multi-arterial LAVA acquisitions are given in (b,c,d). Last, LAVA-Flex is acquired again during the portal phase (e)

Figure 3 . 8 :

 38 Figure 3.8: Placement of the ROIs on the three hepatic veins (right, middle and left) on the three LAVA images acquired during the arterial phase.

Figure 3 . 9 :

 39 Figure 3.9: The values of hepatic veins enhancement at the arterial time are plotted as a function of the lung shunt percentage for every patient whose CT-scan was analysed. The linear regression line that best fits the data is also shown. The correlation coefficient value is 0.68±0.02. Except for some patients that appear to be unrelated, the graph shows a correlation between the two variables

Figure 3 . 10 :

 310 Figure 3.10: The ratio between the hepatic veins enhancements at the portal time (∆C p ) is plotted as a function of on the relative lung shunt fraction. The correlation coefficient value is 0.69±0.01

Figure 3 .

 3 Figure 3.11: The difference in signal intensity between the two ROIs placed on the vena cava before and after the arrival of hepatic veins is plotted based on the relative lung shunt fraction. The correlation coefficient value is 0.29±0.16

( a )

 a Figure 3.12: Enhancement curves of the three hepatic veins (left LHV, median MHV and right RHV) during the arterial phase, realised through the analysis of the multi-arterial MRI acquisitions.

Figure 3 .

 3 Figure 3.13: CT image acquired during the arterial phase. The presence of an extensive tumour prevents the localisation of the mean hepatic vein. Other hepatic veins are marked by yellow ROIs

Fig. 3 .

 3 14 shows the comparison between the first LAVA sequence, with a better signal to noise ratio but a too long acquisition time, and the rapid one tested on the same volunteer.

( a )

 a Figure 3.14: Comparison between the two LAVA sequence tested on a healthy volunteer. The first image has a longer acquisition time but a very good signal to noise ratio. Hepatic veins are more visible due to the better spatial resolution. The right image shows the rapid LAVA sequence, allowing to capture an image every 9 seconds, that can be used for multi-arterial coronal imaging

Figure 4 . 1 :

 41 Figure 4.1: Outline of the method used for a patient-specific simulation of the hepatic arterial tree. In order to initialise the model, segmentation of liver with tumours (from portal-phase images) and arteries (from angiography) is needed

Fig. 4 .

 4 Fig. 4.2a illustrates a liver MRI of a patient suffering from cirrhosis and HCC, acquired during the portal phase, offering the maximum parenchyma enhancement. It can be seen that, in addition to an irregular shape due to cirrhosis, vena cava and portal vein pass through the parenchyma making its borders difficult to distinguish. The presence of a massive tumour complicates as well liver shape and borders. During the portal phase, liver and heart are both enhanced and their similar contrast and close positioning represents another challenge, as can be seen in Fig. 4.2b.

( a )

 a Figure 4.2: Two slices of a MR image acquired at CEM during the portal phase, showing the principal challenges in liver segmentation.

Figure 4 . 3 :

 43 Figure 4.3: Surface evolution during four cycles of the liver shape segmentation in MESA software for Patient 2. Segmentation is based on the active surface method[START_REF] Reska | Fast 3D Segmentation of Hepatic Images Combining Region and Boundary Criteria[END_REF]

Figure 4 . 4 :

 44 Figure 4.4: Example of binary tree ordered in a "depth-first" manner and of its associated input file. Every node can be a bifurcation or a leaf. Leaf nodes are numbered in the order they are met.

Figure 4 . 5 :

 45 Figure 4.5: Example of semi-manual segmentation of vessels with EndoSize. At each step the user defines a new point inside the vessels. Unlike the visual representation proposed in this picture, the vessel centrelines are considered as piecewise straight lines

Figure 4 . 6 :

 46 Figure 4.6: Outline of the workflow leading to the vascular model input file: starting from the 3D angiography, vessels are segmented using a Frangi's filter, before being skeletonised. An analysis of the obtained skeleton finally brings to the definition of a binary tree structure, formalised in the presented manner

  CT-scan at arterial phase (c) CT-scan at portal phase

Figure 4 . 7 :

 47 Figure 4.7: Liver tissue is not enhanced on angiography. Artefacts are due to the presence of metallic coils and surgical clips. White arrow on (b) indicates the visible arteries on arterial phase CT-scan. Arteries are not enhanced on portal phase CT-scan. (Images from CEM)

Figure 4 . 8 :

 48 Figure 4.8: Registration process between portal phase CT-scan and 3D Cone Beam angiography.In order to initialise the model simulating the arterial tree, registration is needed because liver and tumours are segmented from portal phase CT-scan and arteries are segmented from 3D cone-beam angiography

  (a) Part of a binary vascular tree (b) Vascular networks (arteries and veins) and tissue built of macro-cells

Figure 4 . 9 :

 49 Figure 4.9: General overview of the model simulating the hepatic vascular trees and tissue

Figure 4 . 10 :

 410 Figure 4.10: Flow chart representing two loops of events, distinguished in the simulation of the organ (tissue and vascular network) growth

Figure 4 . 11 :

 411 Figure 4.11: Perfusion process of newly created macro-cell, by hepatic vascular network. The optimal configuration of vessels is chosen depending on minimal volume principle. The three vascular trees (hepatic artery, portal vein and hepatic vein) can be simulated by the model. In the context of this thesis, only the arterial tree is considered

  Figure 4.14: 3D Cone-Beam angiography (left) and arteries segmentation in EndoSize (right) for both patients. Non-segmented vessels in (b) are extra-hepatic arteries

  from the angiography of Patient 1 presented all the 46 leaves deduced by the skeleton analysis. In the case of Patient 2, 46 leaves on 78 were found by the algorithm. Two views of the extracted binary tree of Patient 1 are given in Fig. 4.15 (corresponding respectively to the views of Fig. 4.6 and Fig. 4.14a).

Figure 4 . 15 :

 415 Figure 4.15: Example of a binary tree automatically extracted from the arteries segmented with the Frangi's filter on Patient 1

Figure 4 . 16 :

 416 Figure 4.16: Results of liver segmentation in MESA: final 3D surface (upper left) and 2D result (lower right) for 2 patients

•

  Figure 4.17: Registration of arterial CT-scan (right) on portal CT-scan (left). Checkerboard after registration is shown below

Figure 4 . 19 :

 419 Figure 4.19: Registration results between CT-scan and 3D angiography for both patients. Superimposition of liver and tumour contours obtained from portal phase CT and vascular skeleton segmented from angiography

Table 4. 4 :

 4 Input parameters for the simulation of the hepatic arterial tree. Initial number of macro-cells corresponds to the number of leaves of the segmented arterial tree launched on a portable computer equipped with Intel Core i7 CPU (2.3 GHz) and 16 GB of RAM under OS X Yosemite. Computational time for the growth of the entire tree varies between 10 minutes and 1 hour, depending on the number of growth cycles and on liver size. The radii of the simulated vessels are automatically redefined by the model, based on initial pressure and blood flow values. For Patient 1, diameters vary between 2.6 mm at proper hepatic artery and 0.044 mm at the smallest arteriole.Concerning Patient 2, the obtained diameter of the proper hepatic artery is 4 mm, whereas the smallest radius measures 0.06 mm. We recall that the measured diameter of the entrance of the hepatic artery is 4.8 mm for Patient 1 and 6.5 mm for Patient 2. The inconsistency in diameter values will be discussed in the last section if this chapter. Simulated hepatic arterial trees are shown in Fig. 4.20. In Fig. 4.20a only the vessels(a) Only the vessels whose radius is higher than 0.4 mm are shown (b) Complete simulated trees

Figure 4 . 20 :

 420 Figure 4.20: Entire simulated hepatic arterial trees for both patients, including tumour vascularisation

  An optimal catheter tip position is proposed for Patient 1: 71% of the tumour and 3% of healthy tissue are injected (left). Catheter tip is positioned at the root of the arterial tree (right). An optimal catheter tip position is proposed for Patient 2: 64% of the tumour and 2% of healthy tissue are injected (left). Catheter tip is positioned at the root of the arterial tree (right).

Figure 4 . 21 :

 421 Figure 4.21: Results of microspheres distribution simulation, based on Poiseuille's law. Different injection points are compared for the two patients. Light green colour marks arteries where microspheres are present, whereas red arteries are not influenced by the treatment. Blue-yellow colormap shows the percentage of microspheres in the injected tissue. Dark green macro-cells represent tumour tissue where microspheres are absent.

  Figure 4.22: Comparison between simulated and actual microspheres distribution in liver tissue.

  4.19), and the principal simulated ones(Fig. 4.20a). This phenomenon shows the importance of finding a compromise between two constraints: first, to preserve the segmented vessels during the simulation; second, the possibility to modify them during the growth of the arterial tree, in order to verify the optimisation constraints defined in the model when creating new bifurcations. The model has finally been modified in order not to modify the geometrical structure of the input segmented vessels. This new method was applied to the arterial tree shown in Fig. 4.15, automatically extracted from angiography, and results are shown in Fig. 4.23. The image shows that the geometry of the(a) Segmented tree (b) Larger simulated vessels

Figure 4 . 23 :

 423 Figure 4.23: Comparison between the automatically extracted arterial tree of Patient 1 and the larger simulated vessels. The image shows that the geometry of the segmented vessels is generally maintained

  numerical simulations are validated by a physical model of an abdominal aortic aneurysm linked to a Windkessel model, and by phase con-trast MRI on the physical phantom.

  Flow output boundary conditions were determined based on blood flow distribution in hepatic segments estimated from the literature and the model was validated by physical simulations. The 3D blood flow simulation proposed by Aramburu et al. [2016a] used a patient-specific hepatic artery of one patient, made of circular vessels reconstructed from the clinical image. Realistic boundary conditions are computed using a lumped parameter model of peripheral arteries, taking into account standard blood flow distribution in liver segments and flow rate modifications due to the presence of a tumour. A perfusion simplified 2 parameters model is introduced in order to complete a 0-D model version of the flow in the imaged tree. The parameters identification requires the knowledge of synchronised values of pressure and flow rate at the inlet and at different phases of the cardiac cycle (systolic and diastolic phases). The fluid governing equations based on Navier-Stokes equation are here again solved in Ansys Fluent with the FVM. The Lattice Boltzmann method was shown to be a valid alternative to the classical eulerian methods of CFD as finite element and finite volume methods. It easily handles complex geometries of the computational domain and is highly parallelisable. Recent works propose specific 2D and 3D LBM simulations of blood flow in different contexts: among others, Liu [2012] use a biviscosity constitutive relation for the simulation of blood flow in a stenotic artery, Závodszky and Paál [2013] simulate a 3D transient flow of blood in an intracranial aneurysm and Abas et al. [2016a] propose a simulation of 3D laminar flow in an aneurysm geometry comparing three different LB models.

  variables x, p and t denote the three dimensional position in R 3 , momentum in R 3 and time in R. Considering ∆x and ∆p as the finite cubelets centred on x and p in the phase space µ = [(x, p) : x, p ∈ R 3 ], the quantity f (x, p, t)|∆x∆p| gives the expected number of particles in ∆x ∆p at time t, where |∆x∆p| is the Lebesgue measure of ∆x∆p.

Figure 5 . 1 :

 51 Figure 5.1: The two most common 2D and 3D LBM models for the velocity discretisation. Images

Figure 5 . 2 :

 52 Figure 5.2: Examples of extraction of the principal vessels of the simulated tree, according to different diameter criteria

Figure 5 . 3 :

 53 Figure 5.3: LBM blood simulation in the patient-specific arterial tree of Patient 1. Vessels with diameter > 2 mm represent the simulation domain. On the left, the vessels borders, where velocity is zero due to no-slip conditions, are visible. On the right image, regular slices of the vessels illustrate velocity at the interior of the vessels

Figure 5 . 4 :

 54 Figure 5.4: LBM blood simulation in the patient-specific arterial tree of Patient 2. Vessels with diameter > 2 mm represent the simulation domain. On the left, the vessels borders, where velocity is zero due to no-slip conditions, are visible. On the right image, regular slices of the vessels illustrate velocity at the interior of the vessels

Figure 5 . 5 :

 55 Figure 5.5: An example of the procedure developed for the identification of the input and output surfaces as independent meshes. The vessel ending is first manually cut, and a cap is automatically created in order to distinguish the output surface to the vessels wall in Ansys Fluent

Figures

  Figures 5.6, to 5.9. All the cone-beam angiographies of the 4 patients are shown in Figures 5.6a to 5.9a through a Volume Rendering technique (CEM, Rennes). Patient's tumour and the vessels that irrigate it are clearly visible. It can be noticed in Fig. 5.6 that the vessels supplying

( a )

 a Figure 5.6: Cone-beam angiography and blood flow simulation results on segmented arterial tree for Patient 1

  Figure 5.7: Cone-beam angiography and blood flow simulation results on segmented arterial tree for Patient 2

  Figure 5.9: Cone-beam angiography and blood flow simulation results on segmented arterial tree for Patient 4

Figure 5 . 10 :

 510 Figure 5.10: Comparison between velocity fields obtained with the FVM (left) and the LBM (right), in the extracted hepatic arterial tree of Patient 1. The red arrow indicates the vessel that irrigates the tumour. It can be seen that the velocity is higher in such vessel in both cases. The green arrow indicates the entrance of the arterial tree

Figure 5 .

 5 Figure 5.11: Comparison between velocity fields obtained with the FVM (left) and the LBM (right), in the extracted hepatic arterial tree of Patient 2. The red arrow indicates a restriction of the vessel diameter, probably due to a spasm, where velocity increases. This is not considered in the LBM simulation, due to the different procedure applied for the definition of the computational domain. The green arrow indicates the entrance of the arterial tree

Figure 5 . 12 :

 512 Figure 5.12: Example of complex vascular structures where the parallelised LBM simulation could simulate blood flow. The diameter of the smallest vessels in the presented geometry is 0.5 mm

  Figure 6.2: Scheme of the modified bounce-back rule defining the update of the density distribution functions around the i boundary link at the interface between solid and fluid nodes. Image taken

Figure 6 . 3 :

 63 Figure 6.3: LBM blood simulation in a 2D representative arterial geometry

Figure 6 . 5 :

 65 Figure 6.5: Three dimensional LBM simulation of microspheres transport in a 3D representative geometry made of one bifurcation

Figure 6 . 6 :

 66 Figure 6.6: Microspheres trajectories coloured by their velocity magnitude in the artery segmented from the images of the four patients

Figure 6 . 7 :

 67 Figure 6.7: Overview of the principal results obtained for the patient-specific mathematical modelling and simulation of blood flow and microspheres transport

Figure A. 1 :

 1 Figure A.1: At the equilibrium state, the longitudinal magnetisation

Figure A. 2 :

 2 Figure A.2: The signal detected from the MR system is called Free Induction Decay (FID). It corresponds to a sinusoidal wave oscillating at the Larmor frequency, decreasing exponentially according to the constant relaxation time T2 *

(

  FID) and is depicted in Fig. A.2. It corresponds to a sinusoidal wave oscillating at the Larmor frequency, decreasing exponentially according to the relaxation time. During the relaxation phase, two phenomenons occur simultaneously, called T1 relaxation and T2 relaxation. On the one hand, the longitudinal magnetisation → M z grows back towards the

Figure A. 3 :M

 3 Figure A.3: T1 and T2 relaxation curves of grey and white matters. The yellow curve represents the difference between the signal induced by the two tissues: the moment of the acquisition has to be chosen looking for a compromise between a good contrast and a satisfactory signal to noise ratio

  On the other side, the continuity equation (B.11) can be written in a discretised form asN f f J f A f = 0, (B.17) where J f = ρu f is the mass flux through the face f and u f represents the value of u on the cell face f . The pressure-velocity coupling is achieved by introducing a condition linking pressure and velocity, derived from the momentum equation (B.15), into the continuity equation (B.17). In such equation, values of velocity on cells faces, u f , appear. A relationship between face and centre velocity values (u f and u p ) must thus be defined in order to introduce it in the continuity equation. In case of a co-located variable arrangement (as in Ansys Fluent), if velocity on the faces is defined as a linear interpolation of the values on the surrounding cells centres, then an unphysical representation of pressure appears (pressure checkerboard effect). Two solutions to this problem are common in the literature: the use of a staggered grid, where velocity values are stored at cells faces instead of cells centres, or the Rhie and Chow interpolation [Rhie and Chow, 1983], which is generally used in Ansys Fluent. In this method, the a p coefficients of equation (B.15) are used in order to defined a momentumweighted interpolation for face velocity values. Rhie and Chow stated that the discretised momentum equation (B.15) can be also written for the face velocity values. This results in an expression of the face values u f as a function of the coefficient a p , of the velocity values on the centres of the surrounding cells, and of pressure. In Ansys Fluent the flux J f = ρu f

  4. The corrected values of face flux and cell pressure are computed with:p = p * + α p p (B.20)

  Optima MR450w, GE Medical Systems) du Centre Eugène Marquis de Rennes a d'abord été utilisée. La vitesse du fluide étant connue, le calcul de vitesse a été validé. Les contre-indications sont multiples et concernent tout d'abord le dépôt extrahépatique des radiations, notamment dans la vésicule biliaire, les intestins, le duodénum et les poumons. En particulier une irradiation des poumons peut être due à l'existence d'un shunt artério-veineux dans le réseau sanguin hépatique (Fig. 3.1), ce qui permet aux microsphères de rejoindre la veine cave et puis les poumons au lieu de rester bloquées dans les capillaires. Le shunt ne peut pas être traité de manière directe et, dans le cas où les radiations dépassent un certain seuil, il oblige le médecin à exclure le patient du protocole.Le rehaussement des veines hépatiques a été mesuré sur les trois images acquises pendant la phase artérielle, ce qui nous a permis de réaliser, pour chacune des veines, une courbe de rehaussement du contraste. Les résultats affichés sur la figure 3.12 montrent qu'une analyse plus complète d'images IRM comme celles envisagées dans la section 3.7.2, sur une cohorte de patients plus vaste, pourrait mener à une détection particulière du shunt, voire de le quantifier. Le rehaussement d'une veine hépatique en particulier, comme pour le cas de la figure 3.12a, montre qu'une courbe de rehaussement des trois veines peut apporter des informations beaucoup plus précises et approfondies qu'une seule acquisition TDM. Une autre méthode, plus automatique, a été développée ensuite: un filtre de Frangi est appliqué à l'image d'angiographie dans le but de rehausser les vaisseaux, qui sont finalement segmentés suite à un seuillage manuel et une division automatique en composantes connexes. L'arbre obtenu est soumis à une squelettisation et ensuite analysé de façon automatique dans le but de définir une structure d'arbre binaire.Le foie et les lésions tumorales ont été segmentés avec le logiciel Mesa, qui se base sur un modèle de surface active. Une étape de recalage se révèle nécessaire dans la phase d'initialisation du modèle de croissance vasculaire, afin de placer les vaisseaux segmentés, ainsi que la forme du foie, dans un même repère. Cela exploite la TDM à la phase artérielle des vaisseaux segmentés par le filtre de Frangi (cf. Chapitre 4). Les images binaires obtenues grâce à la segmentation ont été utilisées pour construire de façon semi-automatique un maillage correspondant à la surface des vaisseaux et aux surfaces d'entrée et de sortie. Les mêmes hypothèses que dans la méthode précédente ont été posées ici sur les caractéristiques du sang. La vitesse initiale à l'entrée de l'arbre artériel a été définie à partir de données de la littérature[START_REF] Hubner | Hepatic arterial blood flow velocities: assessment by transcutaneous and intravascular Doppler sonography[END_REF]. Étant donné la nature rétrospective de cette étude, il n'a pas été possible de mesurer la vitesse par contraste de phase sur les mêmes patients. En sortie de l'arbre, le débit a été défini, de façon simplifiée, comme étant proportionnel au diamètre de chaque vaisseau. Ce choix a été conforté par Pour ce qui concerne la densité de la suspension (diluée ou concentrée), cela dépend évidemment du rayon du vaisseau. Le diamètre des microsphères injectées varie entre 20 et 30 µm, et leur nombre peut arriver jusqu'à 8 millions de sphères. On peut donc calculer que la distance moyenne entre une microsphère et une autre, du moins pendant l'injection, correspond à environ cinq fois la taille d'une microsphère. Une fois injectées dans le sang, la fraction volumique devient encore plus petite (de l'ordre de 5 10 -5 ), et le diamètre des microsphères est plus que 100 fois plus petit que le diamètre du vaisseau. Cela permet de considérer la suspension comme étant diluée. Dans les plus grands vaisseaux donc, l'interaction entre sphères ainsi que l'influence des sphères sur le sang peut être négligée.

		Les deux séquences optimisées ont finalement été introduites dans le protocole clin-
		ique et réalisées jusqu'à présent sur 8 patients atteints d'un cancer du foie. Les mesures En conclusion, un modèle d'arbre artériel hépatique spécifique-patient, permettant une
	La séquence mesure la vitesse des protons dans les trois directions spatiales et propose une sont en général possibles seulement dans les premières bifurcations de l'artère hépatique, simulation des vaisseaux jusqu'à un diamètre d'environ 0.02 mm, a été développé. Les La faible disponibilité de données réelles pour définir les conditions aux limites est un
	image affichant la norme des trois composantes. jusqu'à environ 2 mm de diamètre. L'anatomie de l'artère hépatique ayant une très grande valeurs de pression et de débit utilisées dans l'initialisation du modèle devront être adap-problème majeur dans la littérature de la mécanique des fluides. Une approche envisagée
	Des acquisitions sur six volontaires sains ont ensuite été réalisées dans le but d'optimiser variabilité inter-patients, il n'a pas toujours été possible d'identifier la bifurcation entre tées aux particularités du patient, notamment par des mesures de pression dans l'arbre consiste à coupler le modèle simplifié d'arbre entier, avec la simulation par CFD, plus
	les paramètres de la séquence. Le tuyau de l'objet test a aussi été placé sous le dos du artère hépatique gauche et droite, voir même l'artère hépatique propre. Lorsque cela a été artériel, qui sont envisagées. Le modèle devra être amélioré dans le but de conserver plus réaliste, dans l'arbre initial segmenté. Dans ce contexte, les valeurs de pression et débit
	volontaire en tant que référence. possible, les vitesses de l'aorte abdominale, la veine cave, la veine porte et l'artère rénale d'informations de l'arbre initial, notamment le diamètre et la géométrie initiaux. disponible partout dans l'arbre complet permettraient d'initialiser avec des conditions aux
	Des séquences 3D ont d'abord été testées dans le but d'imager le plus possible de ont été mesurées en tant que référence. Les mesures de vitesse moyenne sont affichées dans Une identification précoce et non invasive de cette contre-indication permettrait au patient limites plus réalistes la simulation numérique par CFD.
	vaisseaux. Une synchronisation cardiaque de séquences 3D (3D + t, aussi appelée 4D) le Tableau 2.5. Des images d'IRM PC en 3D acquises sur patient sont présentées sur la d'être réorienté rapidement vers d'autres solutions thérapeutiques et de pouvoir améliorer 5 Modélisation spécifique-patient du flux sanguin
	Fig. 2.19 et Fig. 2.20. son espérance de vie.
		Les images cliniques des patients ont d'abord été analysées dans le but d'extraire des
		données pour l'initialisation du modèle. L'angiographie rotationnelle 3D, acquise pendant
		l'artériographie diagnostique, offre un portrait très détaillé de l'arbre vasculaire injecté avec
	6	le produit de contraste, jusqu'à un diamètre d'environ 0.5 mm. Le parenchyme hépatique Modélisation spécifique-patient du transport des microsphères
		est rehaussé pendant la phase portale de la tomodensitométrie, le foie et les tumeurs sont Les microsphères radioactives sont injectées dans l'artère hépatique diluées dans au moins
	donc segmentés à partir de cette acquisition. 20 cc de sérum physiologique. Étant donné la taille des microsphères, la suspension n'a

2.11). Ainsi, le contraste de l'image est déterminé par le déphasage des spins, qui est proportionnel à leur vitesse. Pendant le premier lobe du premier gradient bipolaire la phase des spins augmente, pendant le deuxième lobe le déphasage est annulé, sauf pour les protons en mouvement, où les déphasages s'additionnent (Fig.

2

.10). Le deuxième gradient est appliqué avec une polarité inversée par rapport au premier. Les deux images sont alors soustraites dans le but de réduire le signal dérivant des tissus immobiles et de doubler celui des protons en mouvement. Le déphasage maximum étant de ±180°, une vitesse d'encodage maximale existe et doit être fixée comme paramètre de l'acquisition en fonction de la vitesse du flux étudié, idéalement comme étant légèrement supérieure.

L'imagerie par contraste de phase est sujette à des artefacts importants, surtout dans la direction d'encodage de la phase. Ceux-ci peuvent être limités en adaptant le choix des paramètres d'acquisition comme la vitesse d'encodage, le champ de vue et la direction d'encodage de la phase, dans le but d'optimiser le rapport signal sur bruit.

Des acquisitions sur objet test ont d'abord été réalisées dans le but de valider le calcul de vitesse à partir du niveau de gris de l'image acquise. Un tuyau en plastique de diamètre interne de 3 mm a été connecté à l'injecteur de produit de contraste et du sérum physiologique a été injecté avec un débit de 0.4 ml/s. L'objet test a été placé de façon à pouvoir mesurer le flux dans les trois directions. Une séquence 3D PC disponible sur la machine IRM 1.5 T (n'étant pas disponible, seulement des mesures de vitesse moyenne ont été obtenues.

Différentes valeurs de vitesse d'encodage ont été testées, et une valeur de 50 cm/s s'est révélée adaptée pour recouvrir le spectre de vitesses dans l'artère hépatique, en considérant la vitesse plus importante pendant la phase systolique. Dans le but d'améliorer le rapport signal sur bruit mais de réduire les effets de volume partiels, une épaisseur de coupe de 2.4 mm a été finalement sélectionnée, ainsi qu'une taille de pixel de 0.94 mm. Étant donné la très grande variabilité de la direction des vaisseaux, un angle de bascule θ et un temps de répétition TR relativement petits ont dû être choisis et les valeurs de θ = 10 et TR = 10.6 ms ont été retenues. Une valeur de 20.8 kHz pour la bande passante représentait un bon compromis entre la durée de l'acquisition et le rapport signal sur bruit. À cause de la géométrie de l'artère hépatique (Fig.

1

.4), une acquisition dans le plan coronal offrait une meilleure visualisation du vaisseau. La direction d'encodage de la phase produisant le moins d'artefacts était la direction droite-gauche. Le nombre de coupes a ainsi pu être réduit à 30 dans le but de limiter la durée de l'acquisition, sans pour autant induire des artefacts importants dans la région d'entrée du foie.

Nous avons décidé de réaliser aussi des séquences en 2D, dont la durée permet une acquisition en apnée, et donc avec moins d'artefacts liés au mouvement. Le plan sagittal, en général perpendiculaire à la direction de l'artère hépatique, rejoignant le foie à partir de l'aorte, a été retenu. Cela offre plus de liberté aux manipulateurs dans la sélection du plan, en augmentant les chances d'imager ce vaisseau sans effort particulier. Le choix d'un plan perpendiculaire au vaisseau a aussi permis d'augmenter l'épaisseur de coupe à 8 mm et l'angle de bascule à 70°, et d'obtenir ainsi un meilleur rapport signal sur bruit. Par contre, l'absence de coupes adjacentes complique le repérage des vaisseaux, qui peut se révéler très difficile, même suite à une comparaison avec d'autres images de la même région. La validation de telles mesures n'est pas une tâche évidente étant donné le faible nombre de techniques disponibles pour obtenir une vérité terrain. La présence de l'objet test placé sous le dos du patient nous a permis de valider la technique d'acquisition. Une échographie Identification des patients à risque Depuis le début du traitement par Theraspheres ® au Centre Eugène Marquis, environ la moitié des patients qui ont été inclus dans le protocole n'a pas pu finalement bénéficier du traitement. Les données des patients qui ont été inclus dans le protocole des Theraspheres ® au CEM dans les dernières années ont été recueillies et analysées dans le but de déterminer des paramètres indicateurs d'une présence de shunt, mesurables de façon précoce et non invasive sur les images TDM ou IRM du patient. L'hypothèse d'une corrélation entre le pourcentage de shunt pulmonaire et le rehaussement précoce des veines hépatique après injection du produit de contraste a été posée. Le rehaussement du contraste a ainsi été mesuré dans les veines hépatiques sur les images acquises avant l'injection, pendant la phase artérielle et la phase portale. Dans le but d'analyser la cinétique du produit de contraste, notamment l'origine d'un possible rehaussement précoce dans la veine cave, deux mesures ont aussi été effectuées en aval et en amont des veines hépatiques, dans la veine cave, pendant la phase artérielle (Fig. 3.6). Une analyse statistique a été réalisée sur les données extraites à partir des images TDM de 59 patients. Une corrélation a été mise en évidence en particulier entre le pourcentage de shunt et le rehaussement précoce des veines hépatiques (Fig. 3.9). Un test de Student a été appliqué à un échantillon comprenant les 29 patients avec shunt et les patients sans shunt, donnant lieu à des résultats prometteurs notamment concernant le rehaussement des veines hépatiques (p-value < 0.05). Les images IRM offrent une bien meilleure visibilité du réseau sanguin, ainsi qu'une plus fine résolution temporelle que le scanner. La possibilité de réaliser de l'imagerie rapide en IRM dynamique permet un suivi détaillé de la propagation du produit de contraste. Malgré le faible nombre de patients avec shunt ayant subi un examen par IRM, leurs images ont été analysées en suivant les mêmes critères. Sur 9 des 18 patients il a été possible d'identifier et mesurer les trois veines hépatiques, en donnant lieu à une analyse distincte pour chacune. 4 Modélisation spécifique-patient de l'arbre artériel hépatique Un modèle de croissance vasculaire précédemment développé au sein de l'équipe a été adapté dans le but de réaliser une simulation spécifique à chaque patient. À partir de quelques vaisseaux initiaux, de la forme du foie et éventuellement de régions tumorales, le modèle simule un arbre complet sur la base de conditions de pression et/ou débit en entrée et en sortie de l'arbre. La naissance de nouvelles bifurcations de l'arbre est régulée par la loi de conservation de la matière et par une relation algébrique entre le rayon du vaisseau père et des vaisseaux fils. Dans chaque branche de vaisseau comprise entre deux bifurcations, la loi de Poiseuille est utilisée pour simuler de façon simplifiée l'écoulement du sang. Cette relation permet de relier la géométrie du vaisseau au débit et à la pression du sang. Des caractéristiques différentes concernant la vascularisation et le débit sanguin sont affectées aux régions tumorales. Le tissu hépatique est considéré comme un ensemble de macro-cellules virtuelles, d'une taille d'environ 0.125 cm 3 , auxquelles sont affectées certaines propriétés de débit et de densité, ainsi que de probabilité de naissance et de mort. Cet ensemble de macro-cellules est alimenté par l'arbre vasculaire simulé. Dans le modèle vasculaire, l'artère hépatique est considérée comme un arbre binaire, et une telle structure doit donc être reconstruite à partir de l'image d'angiographie. Dans un premier temps, cela a été fait de façon semi-automatique avec le logiciel Endosize © (Therenva SAS, Rennes, France). pour effectuer une étape de recalage intermédiaire. En effet, il s'agit de la seule image où la forme du foie et les plus grands vaisseaux sont visibles. Un recalage multimodal automatique est utilisé pour recaler de façon rigide la TDM artérielle sur la portale. Une segmentation grossière des plus grands vaisseaux est utilisée pour recaler l'angiographie 3D sur la TDM artérielle. Une méthode de recalage multimodale est ensuite utilisée pour identifier de façon automatique la matrice finale de transformation affine, permettant le recalage des vaisseaux issus de l'angiographie avec la forme du foie issue de la TDM portale. Une fois le modèle vasculaire correctement initialisé avec les données géométriques du patient, un arbre adulte spécifique-patient a été obtenu pour deux patients. Les valeurs de pression en entrée et en sortie ont été estimées grâce aux données de la littérature. Les arbres obtenus semblent bien représenter, au moins d'un point de vue statistique, les arbres des patients et sont visualisés sur la Fig. 4.20. A ce niveau, une première simulation, simplifiée, de l'injection des microsphères a été introduite dans le modèle: il est supposé ici que les microsphères se répartissent de façon uniforme en aval du point d'injection, proportionnellement aux débits donnés par la loi de Poiseuille. La distribution des microsphères peut alors être analysée pour différents points d'injection, correspondant à différentes bifurcations dans l'arbre. Il est alors possible de quantifier le pourcentage de tumeur et de tissu sain injecté (c'est-à-dire où des microsphères sont présentes), ainsi que le pourcentage de microsphères qui rejoignent la tumeur et celles qui vont dans le tissu sain. Une vérité terrain qui décrit de façon précise la distribution des microsphères étant possible seulement par histologie, la scintigraphie post traitement est utilisée pour estimer leur distribution. Une comparaison entre le pourcentage d'activité déposée dans la tumeur au cours de la simulation et celui mesuré en clinique montre une prévision raisonnable de notre modèle (Tableau 4.5). Dans le Chapitre 5, deux méthodes pour la simulation du sang dans un arbre spécifique patient sont comparées: la méthode de lattice Boltzmann et la méthode des volumes finis. Ces deux méthodes ont largement été utilisées dans la littérature pour des simulations de l'écoulement sanguin. Des simulations spécifique-patient dans l'arbre artériel hépatique sont beaucoup plus rares, étant donnée la complexité des vaisseaux, leur petite taille et la difficulté d'obtenir des mesures fiables du flux sanguin, et des pressions. La méthode de lattice Boltzmann est connue pour sa capacité à manipuler des domaines avec des géométries complexes et pour sa parallélisation aisée. Cette technique dérive de la théorie cinétique des gaz: le fluide est considéré comme un ensemble de particules mésoscopique qui se déplacent sur un réseau prédéfini. Les particules sont contraintes à bouger selon un ensemble discret de vitesses et leur collision est gouvernée par la conservation de la masse et du moment. L'évolution du système est décrite par l'équation de Boltzmann discrétisée. Cela décrit la propagation et la collision des particules par une fonction de distribution f i (x, t) qui représente la probabilité de trouver au temps t une particule au noeud x ayant une vitesse e i . La méthode, précédemment développée pour la simulation dans un arbre représentatif, a été utilisée pour la simulation du sang dans les deux géométries spécifique-patient générée dans le chapitre précédant. Le sang a été considéré de façon simplifié comme étant un fluide newtonien, avec une densité de 1051 kg/m 3 et une viscosité cinématique de 3.8 10 -6 m 2 /s. Les conditions aux limites de la simulation ont été définies comme les valeurs de pression données par le modèle de croissance vasculaire (cf. Chapitre 4). La carte de vitesses obtenue dans l'arbre vasculaire d'un patient est donnée sur la Fig. 5.3. En ce qui concerne la méthode des volumes finis, le logiciel Ansys Fluent a été utilisé pour résoudre les équations de Navier-Stokes. Le domaine de simulation a ici été extrait directement à partir une analyse du rapport entre débits et diamètres des vaisseaux, réalisée sur un volontaire à partir des images IRM PC. Les équations ont été résolues avec un algorithme "SIMPLE". Les résultats sont présentés sur les figures 5.6 à 5.9. Les résultats obtenus dans les premières bifurcations de l'arbre hépatique ont été validés de façon qualitative par comparaison avec les valeurs de vitesse mesurés par IRM PC sur d'autres patients, et par les données de la littérature. Les valeurs de vitesse et pression plus loin dans l'arbre hépatique sont très peu connues; seules des valeurs génériques concernant tout le circuit systémique sont disponibles et semblent être cohérentes avec nos résultats. Une étude clinique complète devra être effectuée, dans le but de réaliser toutes les acquisitions nécessaires pour un même patient. Un réglage des conditions aux bords pourra donc être réalisé, dans le but de les adapter aux valeurs de vitesse mesurées à differents endroits de l'arbre. De plus, une acquisition PC 3D + t sera réalisée dès que ce type d'acquisition sera disponible sur l'imageur IRM du Centre Eugène Marquis, afin de tenir compte des différentes phases du cycle cardiaque. La segmentation des vaisseaux devra aussi être adaptée aux caractéristiques de chaque image (présence de bruit, complexité des vaisseaux irriguant la tumeur, ...) pour obtenir une représentation plus précise du profil des vaisseaux. pas été considérée comme une deuxième phase, solide, du flux, mais une description plus précise des sphères a été considérée, de façon à pouvoir explicitement calculer la force de trainée exercée sur chaque sphère par le sang. La simulation des microsphères avec la méthode LBM a été développée de façon à pouvoir considérer le transport des sphères même dans les petits vaisseaux, là où l'influence des sphères sur le fluide doit être prise en compte. Dans cette méthode, les microsphères sont considérées comme un ensemble de noeuds du réseau (le diamètre chaque micro-sphère est constitué d'au moins dix noeuds). Dans les noeuds à l'interface entre le fluide et les sphères, l'évolution des fonctions de distribution f i (x, t) des particules entrant en collision est gouvernée par le principe de conservation de la quantité de mouvement. L'échange de quantité de mouvement à l'interface entre le fluide et une sphère permet aussi de calculer la force hydrodynamique exercée par le fluide sur les sphères. La collision des sphères entre elles et avec la paroi du vaisseau est simulée de manière parfaitement élastique. L'équation du principe fondamental de la dynamique est donc résolue pour chaque microsphère. Une méthode des différences finies centrées est enfin utilisé pour calculer la vitesse des microsphères. La méthode a été pour l'instant appliquée sur une géométrie simple at avec une configuration simplifiée des sphères; les résultats sont présentés sur les Fig. 6.4 et 6.5. Le transport des microsphères a aussi été considéré par la méthode des volumes finis dans la même géométrie de vaisseaux d'un patient, utilisée dans le chapitre 5 pour la simulation du flux sanguin. Dans cette géométrie, qui ne représente que les plus grands vaisseaux de l'arbre artériel hépatique, le rapport entre le diamètre des microsphères et le diamètre du vaisseau est très faible, et l'interaction entre sphères ainsi que et l'influence des sphères sur le fluide ont donc été négligées. Une simulation de type "one-way" (influence du fluide sur les sphères uniquement) a donc été utilisée dans le logiciel Fluent, ce qui a permis de calculer les trajectoires des microsphères dans un temps de calcul négligeable, une fois le champ de vitesse du sang calculé. Les résultats sont affichés sur la Fig. 6.6. Ce travail ouvre la voie vers une modélisation spécifique-patient de la RTMI. La plupart des étapes, de l'analyse du protocole d'imagerie jusqu'à la simulation des trajectoires des microsphères dans l'arbre artériel, ont été étudiées à travers une approche basée image. Des méthodes basées sur l'IRM et l'extractions de données fonctionnelles du patient ont aussi été développées. Dans le but de compléter le caractère spécifique-patient du modèle, la priorité des travaux futurs, portera sur l'initialisation des modèles de CFD avec les données de flux spécifique patient.

  Venc_scale being the velocity encode scale, given in the dicom header of the image. As shown in Fig.2.14, in the case of the Optima MR450w GE system of CEM, the velocity2.2. PROPOSED IMPROVEMENTS OF THE CURRENT CLINICAL PROTOCOL 81Table2.1: Parameters of the 3D phase contrast sequence used to measure fluid velocity on the

	test object
	Acquisition type	3D
	Velocity encoding 30 cm/s
	Slice thickness	2 mm
	Number of slices	160
	Pixel size	0.43 mm × 0.43 mm
	Field of view	220 mm × 220 mm
	Repetition time	9.9 ms
	Echo time	4.8 ms
	Flip angle	15°A
	cquisition time	4:02 min
	value / Magnitude pixel value Vscale	, Vscale = π	Venc_scale Venc	,	(2.1)

Table 2 . 2 :

 22 Parameters of the optimised 2D and 3D phase contrast sequences we propose for the quantification of blood velocity in the hepatic artery

	Acquisition type	2D	3D
	Velocity encoding	50 cm/s	50 cm/s
	Slice thickness	8 mm	2.4 mm
	Number of slices	1	60 (interpolated from 30 acquisitions)
	Pixel size	0.94 mm × 0.94 mm 0.90 mm × 0.90 mm
	Field of view	240 mm × 240 mm	460 mm × 460 mm × 72 mm
	Repetition time	17 ms	10.6 ms
	Echo time	6.7 ms	4.5 ms
	Flip angle	70°10°R
	eceiver bandwidth	15.6 kHz	20.8 kHz
	Acquisition plane	sagittal	coronal
	Phase encoding direction heat-to-feet	right-to-left
	Acquisition time	18 s	1:23 min

Table 2 .

 2 3: Blood velocity measured on the proposed optimised 3D PC MRI sequence on a healthy

	volunteer
	Aorta	18 cm/s
	Vena cava	19 cm/s
	Celiac trunk	27 cm/s
	Common hepatic artery 5 cm/s
	Left hepatic artery	6 cm/s

Table 2 .

 2 4: Blood velocity measured on the proposed optimised 2D PC MRI sequence on a healthy

	volunteer
	Renal artery	24 cm/s
	Portal vein	15 cm/s
	Celiac trunk	29 cm/s
	Hepatic artery 18 cm/s
	Test object	7 cm/s

Table 2 .

 2 5: Velocity measurements realised on the 3D and 2D PC optimised sequenced acquired on 8 HCC patients (P1, ..., P8) at CEM and the average value. HA stands for hepatic artery, PV for portal vein and VC for vena cava. Values are given in cm/s

		Left HA Right HA Proper HA Common HA Celiac		PV	Renal artery Aorta VC
							trunk
		3D	3D	3D	3D	2D	3D	3D 2D 3D	2D	3D	3D
	P1	11	7	13	15	19	37	15 20 35	36	32	14
	P2	-	-	-	12	26	49	11 12 35	-	26	30
	P3	13	17	22	23	36	25	6	8	-	-	-	20
	P4	-	-	12	23	43	25	17 25	-	-	-	-
	P5	8	14	25	22	29	26	15 14	-	-	-	21
	P6	12	13	16	24	28	37	11 13 23	30	-	14
	P7	13	8	15	23	38	31	11 21	-	-	28	-
	P8	7	5	8	9	20	16	8	13 29	30	13	15
	Avg	11	11	16	19	30	31	12 16 30	32	25	19

Table 2 .

 2 6: Comparison of average blood velocities in aorta and vena cava measured with US, 2D PC MRI and 3D PC MRI on the same volunteer

		Ultrasounds 3D PC MRI 2D PC MRI
	Aorta	20 cm/s	16 cm/s	22 cm/s
	Vena cava	15 cm/s	25 cm/s	29 cm/s

Table 3 . 2 :

 32 The table gives the correlation coefficient and the parameters of linear regression approximation concerning all the values measured on the CT-scan data

	Considered value Correlation coefficient Equation
	∆C n a	0.68 ± 0.02	∆C n a = 0.03 * Shunt + 0.4
	∆C p	0.69 ± 0.01	∆Cp = -0.01 * Shunt + 0.6
	∆C V C	0.29 ± 0.16	∆C V C = 0.01 * Shunt
	Concerning the statistical tests, F-test was realised in order to decide which Student's
	t-test had to be applied to the different cases. F-test indicated variances inequality in
	the cases of the contrast enhancement of hepatic veins at the arterial phase (∆C n a

Table 3 . 3 :

 33 Results of the F-test and Student's t-test concerning the dependency of the pulmonary shunt percentage on the signal enhancement values concerning the three measured parameters.F-test was used to estimate the equality or not of the variances of the different populations, and consequently to choose which Student's t-test was to be performed. The last column gives the p-values

	measured by the Student's t-test
	Considered value F-test result	Student's t-test p-value
	∆C n a	unequal variances	< 0.001
	∆C p	equal variances	< 0.001
	∆C V C	unequal variances	0.12
	The enhancement values were classified according to the 20% threshold, using the k-

nearest neighbours methods and cross validation, as outlined in Sec. 3.5.2. The obtained table of confusion is shown in Table

3

.4. The presented results lead to a sensitivity value equal to 0.615 and a specificity equal to 0.891. Moreover, the algorithm decision was correct 49/59 times, giving a correct classification percentage equal to 83%. The number of false positives, as well as false negatives is nevertheless important.

Table 3

 3

		Actual class: Shunt > 20% Actual class: Shunt < 20%
	Predicted class: Shunt > 20%	8	5

.4: Confusion matrix obtained for the validation of the classification of the population according to the threshold given by a shunt percentage of 20%. Every patient is classified once: the sum of the values given in the matrix represents the number of patients

Table 3 .

 3 5: Parameters of the rapid multi-arterial LAVA sequence for shunt detection, producing 4 clustered images, every 9 seconds, on the coronal plane

	Acquisition type	3D
	Slice thickness	6 mm
	Number of slices 68
	Pixel size	0.86 mm × 0.86 mm
	Field of view	440 mm × 440 mm × 204 mm
	Repetition time	2.83 ms
	Echo time	1.28 ms
	Flip angle	12°A
	cquisition plane coronal
	Acquisition time 9 s per volume
	compared to the sequence currently employed. Four images were realised every 9 seconds
	during the same acquisition, globally lasting 34 seconds. Currently, one image is separately
	realised just before the injection and three clustered images are obtained next during one
	acquisition: globally, 35.8 seconds elapse for the four images.

Table 4 . 1 :

 41 Geometry of the CT acquisitions used for segmentation of hepatic arteries and liver shape. Values are given in mm

			Field of View	Voxels number	Voxel size
	CT-scan at	Patient 1	407 × 658 × 281	512 × 828 × 281	0.79 × 0.79 × 1
	portal phase	Patient 2	435 × 703 × 335	512 × 827 × 335	0.85 × 0.85 × 1
	CT-scan at	Patient 1	358 × 200 × 274	917 × 512 × 274	0.39 × 0.39 × 1
	arterial phase	Patient 2	577 × 271 × 395	1088 × 512 × 395	0.53 × 0.53 × 1
	Cone Beam	Patient 1	237 × 237 × 172	512 × 512 × 371	0.46 × 0.46 × 0.46
	CT-scan	Patient 2	237 × 237 × 173	512 × 512 × 374	0.46 × 0.46 × 0.46

4.3.1 Image processing results

Table 4 .

 4 2: Frangi's filter parameters for hepatic arteries enhancement on cone-beam angiography

Table 4 . 5 :

 45 Percentage of the injected activity which deposits to the tumour for both patients

		Model computation Ground truth
	Patient 1	9.2 %	9.2 %
	Patient 2	27.8 %	34.7 %

Table 5 .

 5 1: Mesh characteristics for the hepatic arteries of the four patients

		Patient 1 Patient 2 Patient 3 Patient 4
	Number of outlets	24	20	29	50
	Number of faces	2133552	968791	1353653	3995926
	Number of nodes	288175	130792	187330	540732

Table 5

 5

.2 shows area-weighted

Table 5 .

 5

	of the 4 patients
			Patient 1 Patient 2 Patient 3 Patient 4
		Inlet	2.6	2.3	3.6	3.6
	Radius (mm)	Min. of outlets	0.6	0.49	0.5	0.9
		Max. of outlets	2.7	1.3	2.4	2.2
		Average of outlets 1.7 ± 0.6	0.8 ± 0.2	1.2 ± 0.5	1.5 ± 0.3
		Inlet	100.4	106	106	102
	Pressure (mmHg)	Min. of outlets	99.3	96.4	98.9	99. 3
		Max. of outlets	100.3	104.5	101.9	101.7
		Average of outlets 99.8 ± 0.3 101.4 ± 2.1 99.9 ± 0.7 100 ± 0.6
		Inlet	30	30	30	30
	Velocity (cm/s)	Min. of outlets	2.7	11.7	0.03	3.2
		Max. of outlets	4.4	19.8	37	4.8
		Average of outlets 2.8 ± 0.3	12.7 ± 1.7	8.9 ± 6.5	3.3 ± 0.2

2: Radius, area-weighted average pressure and velocity results in segmented arterial trees

Table 5

 5

	.3: Blood flow velocity measured on phase contrast MRI acquisition (pixel size 0.9 mm)
	on a 27 years old healthy volunteer
	Diameter Mean blood Ratio between diameter
			velocity	and flow ratio
	Abdominal aorta	1.5 cm	21 cm/s	-
	Inferior vena cava	1.1 cm	19.4 cm/s	-
	Celiac trunk	7 mm	35 cm/s	-
	Common hepatic artery	3.1 mm	15.5 cm/s	44
	Gastroduodenal artery	3 mm	12 cm/s	59
	Right hepatic artery	2.8 mm	12.9 cm/s	59
	Left hepatic artery	2.2 mm	9.8 cm/s	98
	Splenic artery	3 mm	13.7 cm/s	51

of different sizes (aorta, arteries, arterioles,...) have different blood flow behaviours, this is why only vessels with similar sizes were compared. Even if the velocity measurements on

  Actuellement, leur distribution lors d'une injection est estimée par l'injection préliminaire d'un radiomarqueur, ce qui peut se révéler trop approximatif. Un traitement personnalisé permettrait une concentration des radiations à la tumeur tout en épargnant le tissu sain environnant. Dans ce travail je me suis intéressée au développement d'un modèle numérique, pour une simulation spécifique à chaque patient des trajectoires des microsphères, dans le but d'optimiser le traitement. Le protocole clinique d'imagerie a été exploité et optimisé pour l'extraction de données spécifiques patients telles que le foie, les tumeurs, l'artère hépatique et le flux sanguin. Les tissus et l'artère hépatique (jusqu'à un diamètre de 0.05 mm) sains et malins ont été simulés. Cela nous permet de simuler la distribution des microsphères dans le tissue hépatique, validée grâce à la scintigraphie post-traitement. Il est supposé ici que les microsphères se distribuent de façon proportionnelle au flux sanguin, lequel est modélisé par la loi de Poiseuille. Des simulations plus approfondies en mécanique de fluides numérique du flux sanguin ont ensuite été réalisées dans l'artère hépatique du patient. Pour cela nous avons utilisé et comparé les méthodes des Volumes Finis (Ansys Fluent) et de Lattice Boltzmann (programme développé dans le laboratoire). Le transport des microsphères a été simulé dans l'artère hépatique du patient avec la méthode des volumes finis, et dans une géométrie simplifiée avec la méthode de Lattice Boltzmann. Une séquence IRM de contraste de phase a aussi été optimisée pour l'extraction de la vitesse du sang dans l'artère hépatique, dans le but de valider le modèle numérique.

Doppler a aussi été réalisée sur un volontaire sain. Il a donc été possible de mesurer et comparer la vitesse moyenne du sang dans l'aorte et dans la veine cave, les résultats étant montrés dans le Tableau 2.6.Quelques études proposant des mesures de vitesse dans l'artère hépatique et la veine porte sont disponibles dans la littérature, et sont cohérentes avec notre étude[START_REF] Dyvorne | Quantification of hepatic blood flow using a highresolution phase-contrast MRI sequence with compressed sensing acceleration[END_REF][START_REF] Carlisle | Estimation of total hepatic blood flow by duplex ultrasound[END_REF]. Les mesures réalisées dans l'aorte et la veine cave semblent se différencier plus des valeurs proposées par la littérature et par nos mesures d'échographie. Il est possible que la séquence optimisée pour l'artère hépatique ne soit pas adaptée aux mesures dans des vaisseaux avec des vitesses plus importantes, notamment concernant le choix de vitesse d'encodage.En conclusion, malgré des difficultés concernant la validation des mesures de vitesse et une marge d'erreur importante due à la taille du vaisseau, au mouvement et aux artefacts de l'acquisition, une méthode pour extraire les informations de vitesse à différent endroits de l'arbre artériel hépatique a été proposé. La séquence IRM a été appliquée sur objet test, sur volontaires et a enfin été intégrée dans le protocole clinique, ce qui a une importance capitale pour la réalisation de simulations CFD spécifique-patient.

Remerciements

We thus developed a more efficient procedure. The outline of the new workflow leading to the vascular model input file is sketched in Fig. 4.6.

The result of the segmentation of the vascular structures through a Frangi's filter is first submitted to a 3D skeletonisation process. Even if this phenomenon is quite rare, loops can appear in the skeleton due either to anatomic anomalies or to segmentation.

Loops have been cut according to the lowest intensity voxel criteria. In case of loops due to segmentation, a low voxel intensity is due to a low Frangi's filter output, and thus has the less probability to actually belong to a vessel. Afterwards, an automatic skeleton analysis provides a voxel classification according to three categories: voxels corresponding to branches, to bifurcations and to vessel endings. A "tagged skeleton" is thus obtained and exploited.

We developed an algorithm able to go through the vessels skeleton while extracting a binary tree structure and finally providing the input file of the vascular model. The algorithm first needs the user to define the coordinates of the tree root. Afterwards, it automatically checks the nature of its neighbouring voxels in every direction. The found vessel branch is then followed voxel by voxel. While the algorithm goes through a vessel branch, only a few voxels of the branch are saved in the input file. The spacing between two saved voxels can be manually adjusted based on the vessel curvature. The voxels belonging to a vessel branch are saved on the input file as nodes of class 1 (bifurcation) having only one child. An exemple of this kind of node can be found in the node number 7 of Fig. 4.4. Visited voxels are marked in order to avoid infinite loops. When a bifurcation is found, its data are saved in a temporary structure and its left branch is followed, according to the depth-first order (Fig. 4.4). The algorithm goes through all the left bifurcations until the first leaf is found and its coordinates, class and cell number are saved. The algorithm then jumps to the last visited bifurcation and it follows its right child. The process is repeated until all the left and right children of every bifurcation have been visited.

The automatic skeleton analysis declares as bifurcation every voxel having more than one neighbour. This can lead to some clusters of bifurcation voxels that can be due to segmentation errors and that must be treated with a deeper analysis. Therefore, if it is found that a cluster leads to further tree branches, one of the voxels in the cluster is saved as a bifurcation and the algorithm continues through the found branch. Otherwise, those voxels are neglected.

A temporary input file is finally written. Before the vascular model can be launched, The two values obtained are not necessarily identical. In order to adjust them, the radius of one of the two vessels is modified (c). When one of the two radii has been modified, it is possible to calculate the radius of the father vessel, and then the pressure at the entrance of the bifurcation (d)

at the bifurcation point and in adjusting them by applying a multiplicative coefficient to all the vessels of one of the two subtrees. The radius of the parent vessel is then calculated, and the pressure at the inlet of the structure is deduced.

The growth of a tumour region always corresponds to an increase of the total hepatic DdQq, where d corresponds to the domain dimension and q to the number of the discrete velocities allowed. The most common one in the two dimensional case is the one illustrated in Fig. 5.1a, called D2Q9, meaning that particles are allowed to move according to 9 fixed velocity vectors. In this case the virtual particles can have zero velocity (e 0 ), or move according to the velocities defined below:

e 0 = (0, 0)

(5.10)

In the D2Q9 model the weights w i of the equilibrium distribution function given in Eq.

(5.4) are defined as

(5.11)

The three dimensional model most commonly used, thanks to its stability and computational efficiency, is the D3Q19 model, illustrated in Fig. 5.1b. The 19 allowed velocities are defined as:

e 0 = (0, 0, 0) e 1,2 , e 3,4 , e 5,6 = C (±1, 0, 0) , C (0, ±1, 0) , C (0, 0, ±1) , (5.12) e 7,...,10 , e 11,...14 , e 15,...,18 = C (±1, ±1, 0) , C (±1, 0, ±1) , C (0, ±1, ±1) .

The relative weights of the equilibrium distribution function are:

(5.13)

Boundary conditions

Different kinds of boundary conditions can be introduced in the lattice Boltzmann formulation: velocity, pressure or a combination of both. Periodic, no-slip or classical constant boundary conditions can be implemented, although not explicitly but through the distribution functions f i .

Different methods for the implementation of constant pressure or velocity boundary conditions exist. A very common and simple scheme was proposed by [START_REF] Zou | On pressure and velocity boundary conditions for the lattice Boltzmann BGK model[END_REF] and uses algebraic conditions deduced by the conditions given in the definition of macroscopic All the equations that have been described until now are converted into lattice units for the numerical implementation. In particular, the fluid viscosity in lattice units ν is derived from (5.5) as:

(5.15)

Such equation clarifies even more the strong relationship between the fluid viscosity and the relaxation time, and thus the importance in the choice of the latter. The fluid viscosity in lattice units is also linked to the physical one by:

whereas the lattice density is defined as:

(5.17)

Since ρ 0 is the real fluid density, ρ should be equal to 1. In reality, this is not always verified and ρ is a value close to 1. This is due to the fact that the lattice Boltzmann model does not exactly satisfies the incompressibility condition [START_REF] Feng | Coupled lattice Boltzmann method and discrete element modelling of particle transport in turbulent fluid flows : Computational issues[END_REF].

The conversions of space and time units being easily found through a multiplication by h and ∆t, we only specify that the fluid velocity conversion is given by v = C v.

(5.18)

Choice of the implementation parameters

Three fundamental and linked parameters have to be carefully chosen for a stable and precise simulation: ∆t, h and τ . The first two are responsible for the temporal and spatial discretisation error, the last has a strong influence on numerical stability. Some important constraints have to be considered in the choice of such primary parameters.

The most important one is related to a primary hypothesis on which the LBM is based: it has been shown that the solution of the incompressible Navier Stokes equation are recovered from the lattice Boltzmann equation with an error of the order of O(Ma 2 ) [Chen andDoolen, 1998, He andLuo, 1997b]. The method thus provides accurate results only if the computational Mach number verifies Ma 1. In the LBM context, the Mach number is defined as

where v max is the maximum velocity of the simulated fluid and C is the lattice speed. Another important constraint is given by equation (5.5), which says that the three parameters ∆t, h and τ are strongly interconnected and they cannot be chosen independently.

is currently validated by placing a mask on the original image, and whose tuning is still operator-dependent. In particular, the choice of the low threshold needed before the classification into connected components (cf. section 4.1.1) influences the final segmentation especially concerning smaller vessels, whose diameter may be overestimated. In the future an adaptive thresholding method like the Otsu's method [START_REF] Oliveira | Unsupervised retinal vessel segmentation using combined filters[END_REF] will be used in order to extract more realistic vessels diameters and too complex tumour vessels will not be segmented, but they will be simulated by the vascular model described in Chapter 4 instead.

The patient-specific blood flow simulations described in this chapter represent a preliminary, compulsory step that will allow us perform, under the same boundary conditions, a numerical simulation of SIRT microspheres transport in the patient's hepatic arterial tree, thanks to the Fluent particles transport tools and the LB method.

Chapter 6

Patient-specific modelling of microspheres transport

This chapter concerns the analysis and simulation of the suspension of glass microspheres, when they are injected into the hepatic artery and transported by the flowing blood within the hepatic tissues. The trajectories of microspheres are not a priori known and a numerical simulation of this phenomenon could provide to the radiologist a prevision of the distribution of microspheres, in order to optimise the treatment with a better tumour targeting. First of all it is important to understand the phenomenon and to identify the forces that play a role in the transport of Theraspheres ® . Once the physical description of the phenomenon is clarified, the numerical simulation can be addressed. Some of the principal techniques of particles transport simulation will be described hereinafter, before focusing on the ones developed in the context of this work.

The flow of particles in a fluid

The criteria that can be used to classify a two-phase fluid are not absolute and many categorisations exist. [START_REF] Leonardi | Coupled DEM-LBM method for the free-surface simulation of heterogeneous suspensions[END_REF] propose a phenomenological classification of a particulate suspension based on a relative size scale. The classification is not absolute since the characteristics of the phenomena can depend on many factors (particles concentration, ratio between particles size and domain size, fluid viscosity, ...).

For particles considered to have a "small" size, electrostatic interactions are dominant and the particulate suspension follows a Brownian motion. According to Ladd [1994a],

6.2 Lattice-Boltzmann method for the simulation of particles transport

The lattice Boltzmann method has been improved over the years for the simulation of multi-component fluids. The simplest version consists in an improvement of the classical Bhatnagar-Gross-Krook collision operator for the simulation of simple fluids: first, collision between particles of the same fluid are considered, followed by cross collision between the two fluids [START_REF] Asinari | Viscous coupling based lattice Boltzmann model for binary mixtures[END_REF]. Following advances allow the simulation of large mass ratios and ensure thermodynamic consistency [START_REF] Arcidiacono | Simulation of binary mixtures with the lattice Boltzman method[END_REF][START_REF] Arcidiacono | Lattice Boltzmann model for the simulation of multicomponent mixtures[END_REF]. However, simulate a suspension as a two-component fluid means losing precision in the simulation of particles rheology.

More accurate simulation techniques are the point-particle methods and, even more precise, the finite-size particle methods. The first family of methods considers the particle as a mass point with neither size nor shape, and particles interactions are governed by probabilistic laws. For instance, [START_REF] Chopard | Cellular automata and lattice Boltzmann methods: A new approach to computational fluid dynamics and particle transport[END_REF] propose a mixed LBM and Cellular Automata model for the simulation of the transport of snow by the wind. The continuous media is modelled by the Lattice Boltzmann model, whereas the snow transport is simulated by the scattering of point-like particles on a grid, following a law which can be approximated as a Gaussian distribution.

A Lattice Boltzmann method for the simulation of finite-size particles was firstly introduced by Anthony J.C. Ladd in 1994. In his paper [Ladd, 1994a,b], he proposes a new method to integrate into the LBM the description of the interaction between the fluid and moving particles, opening the way to a detailed simulation of every particle's behaviour.

The method, which will be extensively described in section 6.2.1, is based on the conservation of momentum at the interface between the moving particle and the fluid. The hydrodynamic force inducing particles motion is then calculated, based on the momentum exchange. The particle is considered here as a set of lattice nodes, therefore its boundary follows the geometry of the grid, and is approximately considered with a stair-step shape.

Further improvements of such model were made in order to smooth the particle's boundary and consequently reduce instabilities in the computation of the hydrodynamic force.

With this purpose, [START_REF] Noble | A lattice-boltzmann method for partially saturated computational cells[END_REF] proposed an immersed moving boundary technique, further developed by [START_REF] Feng | The immersed boundary-lattice Boltzmann method for solving fluid-particles interaction problems[END_REF], [START_REF] Owen | An efficient framework for fluid -structure interaction using the lattice Boltzmann method and immersed moving boundaries[END_REF] and Leonardi et al. [2012a,b] to cite but a few.

The simulation of turbulent fluid flows was then integrated by [START_REF] Feng | Coupled lattice Boltzmann method and discrete element modelling of particle transport in turbulent fluid flows : Computational issues[END_REF][START_REF] Feng | Combined three-dimensional lattice Boltzmann method and discrete element method for modelling fluid-particle interactions with experimental assessment[END_REF].

The method was further enhanced for the simulation of non Newtonian fluids [START_REF] Leonardi | Numerical rheometry of bulk materials using a power law fluid and the lattice Boltzmann method[END_REF] and of free surface flows [START_REF] Leonardi | Coupled DEM-LBM method for the free-surface simulation of heterogeneous suspensions[END_REF].

The developed method

In order to correctly simulate the injection of glass microspheres in hepatic arteries we decided to perform a simulation which takes into account the flowing of the individual particles, as well as inter-particles and wall-particles collisions.

The decision is justified by numerous reasons. According to [START_REF] Stickel | Fluid Mechanics and Rheology of Dense Suspensions[END_REF],

the simulation of the particles as a continuum is legitimate when the average particles radius is approximately 100 times smaller than the length scale. Obviously, microspheres can be injected even in one of the bigger hepatic arteries, whose diameter can go up to 8 mm [START_REF] Da Silveira | Arterial diameter of the celiac trunk and its branches. Anatomical study[END_REF], but a correct knowledge of the distribution of microspheres in tumour and healthy tissue is provided only if the simulation goes down to small vessels, and our vascular model is able to reproduce vessels down to ∼0.05 mm. In such vessels the ratio between vessel size and particles radii gets largely smaller than 100, being particles radii distribution centred around 0.015 mm. Furthermore, [START_REF] Leonardi | Coupled DEM-LBM method for the free-surface simulation of heterogeneous suspensions[END_REF] suggest that for medium scale particles, a finite-size particle method is more pertinent.

Concerning the density of the suspension of microspheres, this value can be estimated based on the data available on the BTG-Theraspheres ® website [BTG IM -Therasphere

RoW, 2016]. The number of injected microspheres can vary between 1.2 and 8 millions, depending on the needed dose. When released, microspheres are diluted in at least 20 cc of saline solution. Considering the average microspheres diameter as ∼30 µm, we can calculate that the volume ratio between saline and microspheres is around 177, which means that, linearly, the average distance between two adjacent spheres is equal to 5 times the average size of the sphere. This led us to simulate it as a dilute suspension.

We implemented the method firstly proposed by Ladd [1994a,b] and further developed by [START_REF] Feng | Coupled lattice Boltzmann method and discrete element modelling of particle transport in turbulent fluid flows : Computational issues[END_REF] and [START_REF] Leonardi | Coupled DEM-LBM method for the free-surface simulation of heterogeneous suspensions[END_REF] among others. Being such method based on the standard LBM formulation, we will use the same notations than Chapter 5.

Microspheres representation in the lattice First of all, the problem of the definition of a sphere in the lattice is examined. The particle is considered as a set of lattice nodes, as shown in Fig. 6.1. For simplicity's sake, the picture illustrates the 2D case, but the

LATTICE-BOLTZMANN METHOD FOR THE SIMULATION OF PARTICLES TRANSPORT181

same method is applied in the three dimensional context. The shape of the sphere is approximated with a stepwise description. Increasing the number of nodes defining the sphere allows for a more precise depiction. Another consequence of this enhanced precision is the reduction of the fluctuation of boundary shape, during the motion of the particle.

One of the drawback of such fluctuations lies is the development of instabilities during the computation of the hydrodynamic force applied on the interface. Three kinds of nodes are distinguished: the fluid nodes which are adjacent to the spheres, the nodes that belong to the sphere but have at least one fluid neighbour, and the internal nodes of the sphere. Between every internal node and fluid node there is no direct connection. On the contrary, from any of the fluid and solid boundary nodes, there is one of the 8 standard directions (18 in 3D, excluding the centre) which leads to a node of the opposite nature. Such connections are called "boundary links" and their role is fundamental in the description of the interaction between fluid and solid particles.

Interaction between microspheres and fluid The criteria governing the interaction between the transported microspheres and blood is based on the principle of the momentum conservation at the impact between the two components. Specifically, at the interface between the fluid and the particle, the difference between their velocities must be zero. The method describing the momentum exchange between fluid nodes and solid nodes is based on 6.3 The finite volumes method for the simulation of particles transport Different approaches can be used for the simulation of a particles suspension in the finite volumes method. In an Eulerian-Eulerian approach the fluid and the suspension are considered as two interpenetrating continua and the governing equations are solved for both phases. In an Eulerian-Lagrangian approach, the fluid is simulated solving the Navier-Stokes equation, and the particles trajectories are computed by a Lagrangian tracking of every single particle through the flow field.

Based on the size of the particles and their density in the domain, the dominant interaction forces change. When the solid volume fraction is negligible, typically smaller than 10 -6 , and when the interparticle spacing is bigger than 100, the influence of the particles on the flow field is considered to be negligible [START_REF] Utikar | Hydrodynamic Simulation of Cyclone Separators[END_REF]. In classic CFD, this approach is usually defined as "one-way coupling". Under such conditions, the fluid is not modified by the presence of the particles and the particles transport can be computed as a post-processing based on the obtained flow field. A "two-way coupling" is considered when the influence of the particles on the fluid must be taken into account, usually when the volume fraction is bigger than 10 -6 and the interparticle spacing is included between 100 and 10. In this case, the influence of the particles on the fluid should not be neglected.

For denser suspension, the collisions between particles should also be considered. Such conditions are simulated by "four-way coupling" methods.

In an Eulerian-Lagrangian approach, all the forces acting on the particles are taken into account in the Newton second law, that is used to update particles positions and velocities.

The evolution of every single particle is thus described by the following system of coupled ordinary differential equations: If gravity force g is considered in the system, the force F includes the term g(ρ p -ρ) ρ p .

When the values of τ p , u and F are constant over time, the system (6.6) can be solved analytical, giving an explicit expression of u n+1 p and x n+1 p as a function of u n p and x n p . Different numerical schemes can otherwise be used for the computation of particles trajectories. Among others, the Euler implicit scheme and the trapezoidal discretisation scheme. In the latter, the first equation of the system (6.6) is discretised as:

In the trapezoidal scheme, the average values of u * and u * p are defined as:

The particle velocity and position are then updated according to:

The analytic discretisation is very efficient but the trapezoidal one provides an accurate solution in a wider range of conditions.

Method

The microspheres injection was simulated in the same four patient's geometries, created with the procedure described in Chapter 5. The Fluent particles injection tools have been used in order to simulate the trajectories of the radioactive glass microspheres. At the entrance of the hepatic artery, an average diameter of 5 mm can be considered based on our measurements and on the literature. The size of the microspheres is therefore more than 100 times smaller than the artery diameter. In such larger vessels, the microspheres volume fraction is also negligible, of the order of φ 5 10 -5 .

Therefore, a one-way coupling model is used in an Eulerian-Lagrangian approach, neglecting the influence of the microspheres on the fluid and the interparticles collisions.

partition the spatial computational domain in two sections: the principal, segmented tree, down to a certain level, and the smallest vessels simulated with the vascular model. FVM and LBM may be separately used on the two sections: the former on the principal tree, and the latter in the smaller vessels. The coupling between the two simulations would be achieved through an exchange of boundary condition from one section to another. This would lead to optimise the advantages of both methods.

The approach presented here opens the way to wider perspectives in terms of local tumour treatments. Every selective internal therapy, based on the vascular distribution of small vectors, as well as the mock-particles like 99m Tc-MAA, could be optimised through a tailored version of the presented model.

Appendix A

Fundamentals of Magnetic Resonance Imaging

Production of a magnetic field and protons magnetic moment

The operating principles of magnetic resonance imaging are based on the magnetic properties of the atomic nuclei present in the human body. Hydrogen protons have indeed an intrinsic magnetic moment, called spin, due to their continuous spinning around their own axis. Thus, protons act like little magnets, with a North and a South poles. When they are not submitted to a specific magnetic field, their axes are randomly oriented within the tissues. Every proton whose spin is non-zero could be used in magnetic resonance, but due to the high concentration of water and fat in human body, hydrogen protons are used in most applications.

An MR system creates a magnetic field whose strength is considered as an identifier of the machine. 1.5 and 3 teslas are the most common magnetic field strengths of the systems currently used in clinical routine. The three principal MR systems manufacturers (General Electric, Siemens and Philips Healthcare) have proposed 7 teslas systems, but none of them is used for clinical routine yet. When the patient is placed inside the magnetic field, all the hydrogen protons in his body align with the main magnetic field → B 0 induced by the system.

Its axis follows the head-to-foot line and is usually considered the longitudinal direction, also called the z axis. While all spins are aligned to the direction of → B 0 , their sign can vary and some will have the same sign than → B 0 , some the opposite one. As a result, most of the protons will cancel out, but there will be a slight excess of protons aligned in the same sense than → B 0 that will create a net magnetisation, called → M. The presence of the main magnetic z axis following an exponential growth described by:

T1 is the time it takes to this vector to reach the 63% of its original value M z0 , and is a variable associated to every specific tissue and increasing proportionally to the main magnetic field is ∼1000 ms and T2 ∼100 ms. It can be noticed that the moment of the acquisition has to be chosen looking for a compromise between a good tissues contrast and a satisfactory signal to noise ratio.

Definition of an MRI pulse sequence

Spin dephasing inducing the decreasing of the transverse magnetisation → M xy has different reasons. Some of them are random, like their mutual interaction, and some are constant, like magnetic field inhomogeneities. The latter can actually be temporary avoided by the application of a 180°RF pulse. Indeed, if a proton precesses slightly quicker than its neighbour, a 180°RF pulse inverting the phases of the spins will make the faster proton inphase precession, globally allowing a faster acquisition.

Appendix B

The finite volumes method A general conservation law for the scalar quantity u(x, t) can be written under the form:

where F is the flux that determines the transport of u, and f a source term. The flux F can depend on u, on its gradient ∇u or it can be a function of both (second-order equation).

Different schemes can be used for the discretisation of equation (B.1) in time and in space. Time is discretised with a sequence (t n ) n∈N where t 0 = 0 and t n = nk. Thus, a finite difference scheme can be used for the approximation of the time derivative ∂u ∂t (x, t). If a simple explicit Euler scheme is used, this quantity is approximated by ∂u

In order to define a finite volume spatial discretisation scheme, equation (B.1) is integrated over every cell K of the mesh T , giving:

Thanks to the Stokes formula this becomes

where n K (x) is the unit vector normal to the boundary ∂K in x and pointing outwards from K. The second term can be slip into the different parts of the boundary of K:

The quantity F(x, t n ) • n K (x) represents the flux per unit surface area crossing the interface K|L from K to L at x. The flux exchange term

K to L between t n and t n+1 is approximated by a certain quantity F (n) K,L . Using these notations, a finite volume scheme has the form

K is the constant value that approximates the solution u over the time-space cell K × (t n , t n+1 ) .

The conservation law now imposes that on the boundary between every couple of cells K, L ∈ T the flux is conservative:

Defining a finite volume scheme means defining its numerical flux F

K,L as a function of the unknown variable in the surrounding cells. If the flux is a function of the unknown at the time step t n , the method is called explicit. If it is a function of the unknown at the time step t n+1 , then the method is called implicit.

A finite volume method is also characterised by the location, in the mesh, where scalar and vector quantities are stored, whether at the centre of a control volume or on its boundaries. In Ansys Fluent, scalar and vector values are stored by default at the centres of the volumes (co-located variable arrangement), an interpolation of the surrounding values is thus needed for estimating the values on the boundaries that appear in the general equation (B.5).

Discretisation of the Navier-Stokes equations

The Navies Stokes equation are a particular case of the conservation equation (B.2), suitable (and always used in Ansys Fluent) for the diffusion term µ∆u, where no direction has a privileged influence on the transport of the scalar quantity. For the advection term u • ∇u an upwind scheme is more appropriated. Indeed, in such scheme the values u f of velocity components on the cells boundaries are defined as a function of the values of u L at the centres of the cells that are "upstream" to it, relatively to the flux direction. The first-order upwind scheme consists in defining the face value u f as the cell-centre value in the upstream cell: u f = u L . In the second-order upwind scheme, a Taylor expansion of the value u L at the centre of the upstream cell is used: (B.13) where -→ r L is the displacement vector from the upstream cell centre to the cell face centroid.

In order to evaluate the gradients of a scalar quantity, different methods exist. A very efficient one, applied by default in Ansys Fluent, is the least-squares cell-based gradient evaluation method. Defining as r KL the vector from the centre of cell K to the centre of its neighbouring cell L, one can write:

Considering the same equations for all the surrounding cells, the following system is obtained:

where the coefficient matrix J is a function of geometry and ∆u = (u L -u K ), L ∈ N (K).

In order to solve the system, the matrix J is decomposed using a Gram-Schmidt process, yielding to a matrix of weights for each cell. The three components of the weights W x LK , W y LK , W z LK are created for every face K|L, L ∈ N (K) of the cell K. Finally, the gradient ∇u K = ∂u K ∂x , ∂u K ∂y , ∂u K ∂z is computed by multiplying the weight factors W KL by the difference vector (u L -u K ), L ∈ N (K):