Michele Diana

Federico Spadoni

Marco Lorenzi

Christoph Seiler

Marine Breuilly

Alan Garny

Adityo Prakosa

Ján Margeta

Rocio Cabrera Lozoya

Thomas Demarcy

Marc- Michel Rohe

Hugo Talbot

Thomas Benseghir

Pietro Gori

Loic Cadour

Merzieh Kohandani

Mehdi Hadj-Hamou

Nina Miolane

Roch Mollero

Nicolas Duchateau

Aurelie Degletane

Heloise Bleton

Raphael Sivera

Mike Buckingham

Hervé Lombaert

Maxime Sermesant

Xavier Pennec Thank

Frank Blindauer

Mourad Gael Fourré

Chloé Audigier

Invited Talks

Keywords: Modélisation d'ARF, Foie, Personnalisation, Diffusion de la Chaleur, Nécrose Cellulaire, Mécanique des Fluides, Modèle Informatique, Méthode de Lattice Boltzmann, Estimation de Paramètres, Etude Pré-Clinique, Imagerie Medicale RFA Modeling, Liver, Patient-Specific, Heat Transfer, Cellular Necrosis, Computational Fluid Dynamics, Computer Model, Lattice Boltzmann Method, Parameter Estimation, Pre-clinical Study, Medical Imaging v

L'ablation par radiofréquence (ARF) de tumeurs abdominales est rendue difficile par l'absence de planification spécifique à un patient donné. En particulier, l'influence des vaisseaux sanguins et les variations spatiales et temporelles de la conductivité thermique compliquent la prédiction de l'étendue de l'ablation. Cela peut se traduire par des traitements incomplets et l'augmentation des risques de récurrences compromettant le résultat final. En fournissant des outils prédictifs, les modèles biophysiques peuvent aider les cliniciens à planifier et guider efficacement la procédure. Pour que ces modèles soient utiles en routine clinique, le calcul doit être rapide et précis.

Nous présentons d'abord un modèle mathématique détaillé des mécanismes biophysiques impliqués dans l'ARF des tumeurs hépatiques tels que la diffusion de la chaleur, la nécrose cellulaire et le flux sanguin hépatique. Ce modèle simule l'étendue de l'ablation à partir de l'emplacement de la sonde et de certains paramètres biologiques tels que la conductivité ou la capacité calorifique des tissus. Le modèle est basé sur des images médicales, à partir desquelles des modèles personnalisés du foie, des vaisseaux visibles et des tumeurs sont reconstruits. Une nouvelle approche pour résoudre ces équations aux dérivées partielles basée sur la méthode de Lattice Boltzmann est introduite. Le modèle est mis en oeuvre avec une accélération matérielle liée à la carte graphique qui permet d'atteindre un temps de calcul presque égal au temps réel. Les simulations obtenues sont vérifiées précisément avec une solution analytique.

Dans un deuxième temps, nous voulons valider le modèle. Il est d'abord évalué sur des données cliniques de patients qui ont subi une ARF de tumeurs du foie. L'importance de la prise en compte du flux sanguin et de l'estimation des paramètres sont mises en évidence. Ensuite, un protocole expérimental complet combinant des images multi-modales, anatomiques et fonctionnelles, pré et post-opératoires, ainsi que le suivi de la température et de la puissance délivrée pendant l'intervention est présenté. Il permet une validation pré-clinique totale qui considère des données les plus complètes possible pour la validation du modèle.

Enfin, nous estimons automatiquement des paramètres personnalisés en utilisant des algorithmes de résolution de problèmes inverses, afin de mieux prédire l'ablation. Cette stratégie de personnalisation a été testée et validée sur sept ablations dans trois cas cliniques. A partir de l'étude pré-clinique, nous pouvons aller plus loin dans la personnalisation en comparant la température et la puissance délivrée simulées avec la température et la puissance délivrée mesurées pendant la procédure. Cette seconde stratégie de personnalisation a été testée et validée sur douze ablations dans cinq cas pré-cliniques. Ces contributions ont abouti à des résultats prometteurs, et ouvrent de nouvelles perspectives pour planifier et guider l'ARF.

1.1 Motivation

Clinical Context

The World Health Organization (WHO) estimated that 782,000 new liver cancer cases occurred worldwide in 2012, and that it is the second most common cause of death from cancer. The WHO estimated it to be responsible for nearly 745,000 deaths in 2012 that is to say 9.1% of the total cancer death worlwide [Ferlay 2015]. Hepatocellular carcinoma (HCC) for example is one of the most common malignancies encountered throughout the world (more than 1 million cases per year), with increasing frequency in Western countries due to the changing prevalence of hepatitis C [El-Serag 2003]. The prognosis for liver cancer is very poor (overall ratio of mortality to incidence of 0.95) [Ferlay 2015], since treatment of primary and metastatic malignancies in the liver remains a significant challenge in spite of recent advances in cancer therapy. Unfortunately, only less than 25% of patients with primary or secondary liver cancer are candidates for resection or transplantation, which are considered as the most effective treatments. These limitations are due to the patient's condition, the location, size and number of tumors. Proximity to vital organs like vascular and biliary structures is a problem. For example, many patients with cirrhosis or with multiple tumors have inadequate hepatic reserve to tolerate resection [Curley 1997]. Surgical resection are currently restricted to patient with solitary HCC, with liver function extremely well-preserved, without portal hypertension and without abnormal bilirubin. Whilst liver transplantation is limited by the number of donors [Bruix 2002]. Consequently, minimally invasive ablative therapies such as radiofrequency ablation (RFA), high intensity focused ultrasound (HIFU), cryotherapy, microwave tumor coagulation and laser ablation have raised increasing interest for the treatment of liver tumors. Amongst them, RFA is the most widely used approach [Garrean 2008], but it is not yet a reliable alternative in clinical routine. There is a need for training tools for the less-experienced clinicians to improve outcome rates since the success rate increases quickly with experience [Hildebrand 2006]. For the same physician, survival rate of treated patients increased twofolds over a four years period. This learning curve is partly due to the difficult assessment of the cooling effect of the large vessels, porous circulation and blood coagulation, which results in suboptimal ablation and local recurrences in up to 60% of the cases [START_REF] Kim | [END_REF]]. Nevertheless, RFA technology has increasingly been utilized for unresectable hepatic malignancies, is commonly used in the ablation of aberrant conduction pathways in the heart and has become a primary procedure in the treatment of cardiac arrhythmias [Bilchik 2001].

Role of Medical Imaging

Coagulation of living tissue by means of heat was mentioned around 3000 B.C. in the Edwin Smith papyrus, a copy of part of an ancient Egyptian textbook on trauma surgery [Rieder 2013]. It describes 8 cases of tumors or ulcers of the breast that were removed by cauterization with a tool called the fire drill. The writing says about the disease, "There is no treatment"1 . Nowadays, coagulation of living tissue by means of heat is used in clinical routine: radiofrequency ablation (RFA) is an established treatment for liver cancer when resection is not possible. Yet, incomplete treatment and an increased risk of recurrence are common. A tool that would enable the accurate planning of RFA is hence necessary.

To give additional guidance to radiologists or surgeons, several research groups are working on automatic tools to extract information from images and simulate the RFA behaviour in the liver parenchyma. Patient-specific anatomy is not enough since the optimal delivery of RFA is challenged by the presence of large blood vessels and the time-varying thermal conductivity of biological tissue. The effect of the blood flow, the temperature and the cell death on the tissue properties have to be accounted for, thus simulation is necessary. The idea is to combine anatomical data and the knowledge of biological mechanisms in order to build biophysical models that can reproduce the RFA in a patient-specific anatomy and therefore provide predictive tools. By simulating the outputs of the therapy in silico, these models might help clinicians to plan individual treatment for each patient since they would be able to study the effect of the probe position or the heat duration for example on the treatment output.

Theoretical models and computer simulations are powerful tools since they provide vital information on the thermal behavior of ablation rapidly and at low cost. They could provide efficient training tools for the less-experienced clinicians. They could also facilitate the feasibility of new probe geometries, and new RFA protocols.

But, the two most important issues involved in the modeling of RFA include first the experimental validation of the proposed model and secondly the lack of an accurate characterization of the biological tissues [Berjano 2006].

Objectives

Within this scope, the objective of this work is to investigate new computerized predictive tools for patient-specific guidance and therapy planning of radiofrequency ablation (RFA) procedure. Patient-specific simulations may allow clinicians to better understand the intricacy of this procedure and could even allow to rehearse on a virtual patient-specific anatomy. These challenges require many improvements regarding the state of the art. The main questions, that we aim to answer in this thesis, are:

• How to simulate the radiofrequency ablation of hepatic tumors in computational time compatible with clinical settings while taking into account the main occurring biophysical phenomena ?

• Which parameters of this RFA model can be personalized and based on which data ?

• How to validate such RFA model on pre-clinical and clinical data ?

1.2 Main Contributions and Manuscript Organization

Main Contributions

During this PhD study, four main contributions to the modeling of RFA and its validation were proposed.

• A detailed computational model of the biophysical mechanisms (heat transfer, cellular necrosis, hepatic blood flow) involved in RFA of abdominal tumors based on patient images.

• A new implementation of the bio-heat equations coupled with a cellular necrosis model using the Lattice Boltzmann Method on Graphics Processing Units (GPU), which allows near real-time computation.

• The automatic estimation of the main parameters of the model. Two personalization strategies are tested and evaluated on clinical and pre-clinical data.

• The evaluation of the proposed model on clinical data as well as preclinical data from a comprehensive experimental set-up specially designed for RFA model validation.

Organization of the Thesis

The manuscript is organized following the published and submitted works. It follows rather well the studies performed during this thesis in chronological order: from the modeling of RFA to its main parameters personalization and finally to its evaluation on pre-clinical data. This thesis focuses on hepatic tumors and the associated RF ablation procedures. Chapter 2 introduces the background on hepatic anatomy and hepatic functions as well as details on radiofrequency ablation and other minimally invasive ablation techniques. Then, an introduction to the modeling of radiofrequency ablation establishes the state of the art, shows the different models and methods used in the literature, highlight their limits and those further used in our choice of RFA modeling.

Chapter 3 describes our first contribution, published in [Audigier 2013a] and further developed in [Audigier 2015a]. It consists of a new method to compute the extent of ablation based on the Lattice Boltzmann Method (LBM) and patient-specific, pre-operative images. To the best of our knowledge, this is the first time such a complete model, based on medical images has been proposed. A detailed anatomical model of the liver is obtained from volumetric images. Then a computational model of heat diffusion, cellular necrosis, and blood flow through the vessels and liver is employed to compute the extent of ablated tissue given the probe location, ablation duration and biological parameters. Furthermore, a second and improved patient-specific model of blood flow through liver and parenchyma is presented. We rely on the Lattice Boltzmann Method (LBM) to compute not only the heat transfer and cellular death mechanisms but also the blood and parenchyma flow in the liver tissue. This latter method is based on a Computational Fluid Dynamics (CFD) solver which incorporates a porous part to deal with the liver parenchyma. This second framework is particularly efficient for the modeling as it provides a fast solver and naturally accounts for the flow transition between veins and parenchyma. Finally, the model is verified against an analytical solution, showing good consistency.

In Chapter 4, based on [Audigier 2013a] and [Audigier 2015a], we evaluated the predictive power of the proposed framework on ten patients who underwent RFA, for whom pre-and post-operative images were available. Implemented on graphics processing units (GPU), our method simulates 1 minute of ablation in 1.14 minutes, allowing near real-time computation. This represents a 24% speed-up with respect to state of the art computation [Payne 2011]. By providing near real-time computation, our method may enable model-based planning of RFA in clinical settings. Comparisons between the computed ablation extent and ground truth, as observed in post-operative images, were promising (DICE index: 42%, sensitivity: 67%, positive predictive value: 38%). A mean point to mesh error between predicted and actual ablation extent of 10.17 mm is achieved. The clinical criterion is satisfied as the tumor is totally covered by the induced thermal lesion in all cases. The effect of liver perfusion while simulating electrical-heating ablation was also considering. In addition, the results of this chapter demonstrate the need for personalization of tissue properties, as illustrated in one patient. The heat conductivity was adjusted to reduce the point-to-mesh error with the ground-truth from 8.25 mm to 3.85mm for this patient.

Motivated by the results obtained in chapter 4, Chapter 5 presents a framework for parameter estimation in order to model patient-specific radiofrequency ablation of multiple lesions in the case of metastatic diseases, published in [Audigier 2014a]. The proposed forward model is based upon the computational model of heat diffusion, cellular necrosis and blood flow through vessels and liver which relies on patient images presented in chapter 3. We identify the most sensitive material parameters, that need to be personalized from the available clinical imaging and data. The selected parameters are the heat conductivity dt and the porosity ε. These parameters are then estimated using a gradient-free optimization method such that the point-to-mesh distance between the computed necrotic area and observed lesions is minimized. Based on the personalized parameters, the ablation of the remaining lesions are predicted. The framework is applied to a data set of seven lesions from three patients including pre-and post-operative CT images. In each case, the parameters were estimated on one tumor and RFA was simulated on the other tumor(s) using these personalized parameters, assuming the parameters to be spatially invariant within the same patient. Results showed significantly good correlation between predicted and actual ablation extent (average point-to-mesh errors of 4.03 mm).

In order to evaluate the RFA model in a more control environment, a pre-clinical study is proposed in Chapter 6. It is an experimental study especially designed for the validation of subject-specific multi-physics model of liver tumor radiofrequency ablation. The comprehensive experimental set-up is described first. It combines multi-modal, pre-, intra-and post-operative anatomical and functional images, as well as the interventional monitoring of the temperature and delivered power. The RFA model we aim to validate is then introduced. The RFA computation becomes subject-specific after several levels of personalization: anatomic, hemodynamics, heat transfer and an extended cellular necrosis model. Finally, a full processing pipeline required before using the image information into the model is presented. It copes with image noise, resolution and anisotropy. Methods of segmentation, registration and the estimation of smooth vessel tree tailored for computational fluid dynamics simulations are presented. This leads to an end-to-end pre-clinical validation framework that considers the most comprehensive data set for model evaluation.

In Chapter 7, results on the pre-clinical data are presented. The evaluation study encompasses twelve ablations from five healthy pigs: a mean point to mesh error between predicted and actual ablation extent of 5.3 mm is achieved, sufficient for clinical application as it is significantly lower than the difference of 1 cm in diameter configurations of the probe. We evaluate the error induced by each step of the pre-processing. After a sensitivity analysis of the model parameters, we also investigate the estimation of the cellular necrosis and heat transfer models parameters. By fitting the delivered power and the cooling temperatures measured by the RFA probe, key physiological parameters of the model can be correctly estimated as it leads to realistic predictions in terms of necrotic area, power and temperature. Such a framework would therefore not require any necrosis information, and would be better suited for clinical applications.

We conclude on the presented contributions and discuss perspectives in Chapter 8.

Finally, Appendix A gives the derivation of the Lattice Boltzmann equation with the details of the required calculations. The derivations of diffusion equation, advection-diffusion equation and reaction-advection-diffusion equation are studied.

Structure and Function of the Liver

Building an efficient and accurate computational model of RFA requires a high knowledge of the structure and function of the liver. The liver is a vital organ, which has a wide range of functions, including detoxification of various metabolites, protein synthesis, and the production of biochemicals necessary for digestion 2 . A human liver has four lobes of unequal size and shape, whereas pig liver has five lobes. The liver is connected to two large blood vessels, the hepatic artery and the portal vein represented in Figure . 2.1. The hepatic portal vein delivers approximately 75% of the liver's blood supply, and carries blood rich in digested nutrients from the entire gastrointestinal tract, its associated organs and also from the spleen and pancreas. The hepatic artery supply arterial blood to the liver, accounting for the remaining quarter of its blood flow. It carries oxygen-rich blood from the aorta. 3 Oxygen is provided from both sources; approximately half of the liver's oxygen demand is met by the hepatic portal vein, and half is met by the hepatic arteries. Blood flows through the liver sinusoids and empties into each central vein. The central veins coalesce into hepatic veins, which leave the liver and drain into the inferior vena cava. The biliary tree, or biliary tract is derived from the branches of the bile ducts. It is the path by which bile is secreted by the liver then is transported to the first part of the small intestine, the duodenum.

Radiofrequency Ablation

Radiofrequency Ablation Procedure

During RFA procedure, the clinician places one or more probes using percutaneous, laparoscopic or open surgical (celiotomy) techniques within the target area in the liver parenchyma, within the malignant tissue. Although open surgical RFA procedure is more invasive, it allows repeated placement of the RFA probe at different locations to ablate several or large tumors, it can also be performed more easily [Minami 2011]. RF ablation relies on a complete electrical circuit to conduct RF current. This electrical circuit is created from the active probe through the body and is closed by surface electrodes (ground pads) as shown in Figure . 2.2. A complete electrical circuit is then created through the patient body allowing the alternating current of RF energy to go between the needle active electrode and the grounding pads. The white area is the targeted tumor, and the surrounding orange area, the necrosis lesion produced by RFA. Image from [START_REF] Curley 2003steven | Radiofrequency ablation of malignant liver tumors[END_REF]].

From the active electrodes at the tip of the probe, a high frequency alternating electric current flows through the surrounding region and heating is induced due to the electric resistance of the living tissue as illustrated in Figure . 2.3. Depending on the tumor size and the type of RFA probe used, the probe tips are deployed at a certain diameter in order to achieve a desired volume of ablation. This volume is established pre-operatively, ideally to destroy the visible tumor plus a 1 cm safety margin of ablation around it [Lencioni 2005]. Most of the time, in clinical settings, a temperature-controlled mode is used: the delivered RF energy increases up to a given maximum power until the given targeted electrode temperature is reached. Power is then scaled to maintain this temperature during a given ablation time. Following the ablation, the tips of the probe are retracted inside the main axis and the probe tract is cauterized as the RFA needle is withdrawn. This results in thermal coagulative necrosis at temperatures above 50 • C due to irreversible protein denaturation of the cells.

Unlike extirpative therapies, the RFA treatment is difficult to monitor in vivo to prevent recurrence: success of the procedure depends on the complete coverage of the tumor by the thermal zone, which relies on optimal probe placements and the extent of conductive heat delivery controlled by the heating time. However, the latter is challenged by the hepatic blood vessels and the parenchyma perfusion that dissipate heat and make the size and shape of the ablation zone difficult to control thus potentially reducing RFA efficiency and increasing risks of recurrence [START_REF] Kim | [END_REF]]. For example, Figure . 2.4 shows a liver extraction where the effect of the blood vessels on the necrotic area created on an in vivo RFA on a pig liver can be seen. Another great challenge for abdominal tumor treatment with RFA is the pre-interventional planning of the therapy, which has to be accurate in order to achieve a complete destruction of all the tumor cells.

In clinical practice, the planning is based on a 3D CT scan or MR image of the patient acquired pre-operatively. The clinician has to take into account the Besides the poor quality of the U.S guidance, the biophysical aspects have to be accounted for during the intervention. Even if a perfect guidance is achieved, the RFA outcome might still not be optimal. These are the reasons why guidance the RFA protocol for a specific patient is a challenging task and mathematical modeling has the potential to assist the radiofrequency ablation of liver tumors.

Blood vessels

Necrotic lesion

Monitoring of RFA based on MR Thermometry

Currently, the most widespread imaging technique for monitoring RFA of liver tumors is ultrasonography (US). However, real time MR temperature imaging covering the entire tumor and its surrounding area seems to be a more relevant strategy as it could be employed to interactively guide the RF probe to the targeted tumor and to monitor the effect of the therapy [Berjano 2006]. With temperature-based therapies such as RFA, a close relationship exists between cell death and the accumulated thermal dose (TD), which is dependent on both temperature increase and time of exposure [Sapareto 1984]. In their work [Lepetit-Coiffé 2010] show that it is possible to compute voxel-specific TD maps after the suppression of the electromagnetic interference by efficient filtering of the RF signal harmonics at the MR frequency performed like in [Seror 2008], as this otherwise interferes with the measurements of the scanner. They demonstrated in two cases that RF ablation did not result in complete tumor ablation by showing insufficient coverage of the tumors. For ablation under MR guidance, the clinician could be notified in real-time that the end-point of the procedure is not reached whatever the source of the problem (heatsink effect of flow, difficulties in tumor targeting) using MR thermometry. However, despite being the most sensitive imaging modality in detecting liver tumors, MR remains inefficient due to excessive positioning time that results from difficulty in accessing the patient and due to long procedure time [Guenette 2010]. Today, probe thermistors are used to assess the temperature at the tips of the probe. .

Percutaneous Tumor Ablation Techniques

Radiofrequency ablation (RFA) is the most common minimally-invasive technique. However other percutaneous techniques can be used to ablate tumors. We will present most of the currently used techniques in this section.

High Intensity Focused Ultrasound

High Intensity Focused Ultrasound (HIFU) is a medical procedure that applies high-intensity focused ultrasound energy to locally heat and destroy diseased or damaged tissue through ablation 6 . This technique creates short and local elevations of temperature at the focus of an ultrasonic therapy transducer and exposure time on the order of few seconds, inducing irreversible local coagulation necrosis. HIFU has advantages of no-invasiveness, good penetration, good selectiveness and easy power control [Zhang 2006]. Currently, further studies are needed to explore its treatment value for liver or prostate cancer for example.

Microwave Ablation

Microwave Ablation (MWA) is a form of thermal ablation used in interventional radiology to treat cancer. MWA uses electromagnetic waves in the microwave energy spectrum (300 MHz to 300 GHz) to produce tissue-heating effects. The oscillation of polar molecules produces frictional heating, ultimately generating tissue necrosis within solid tumors to damage and kill cancer cells or to make cancer cells more sensitive to the effects of radiation and certain anticancer drugs 7 . Theoretically, MWA has the potential advantages over RFA to generate larger zones of ablation in a shorter time and to generate higher temperatures. However, microwave power is inherently more difficult to generate and deliver safely and efficiently to the tissue when compared to RF power. Currently available microwave systems face technical limitations including large diameter probes, underpowered systems, shaft heating, long and relatively thin ablation zones and unpredictability of the size and shape of the zone of ablation. Further study and continued development of more robust clinical systems is still needed [Lubner 2010].

Cryotherapy

Cryotherapy (also called cryoablation or cryosurgery) is a clinical technique similar to RFA since the energy is delivered directly into the tumor by a probe, which is inserted through the skin. But rather than killing the tumor with heat, cryoablation uses an extremely cold gas (usually argon) to freeze it at a temperature of 233K to 248K. As the gas flows through the needle, a ball of ice crystals forms around the tip of the probe, thus immediately leading to cellular death of the surrounding tissue. This technique has been used for many years by surgeons in the operating room, but in the last few years, the needles have become small enough to be used by interventional radiologists through a small nick in the skin without the need for an operation. This technique has been applied to treat many kinds of tumors, such as breast cancer, primary or metastatic liver neoplasms, renal, lung, pancreas, and prostate cancer [Talbot 2014]. Cryotherapy may be less painful than RFA or MWA and ice ball can be visualized in real time (using CT, MRI or ultrasound) since there is no interference of the procedure with the imaging system. But the set-up is more complicated involving cryogens and tubing, the procedure duration is longer, the maneuverability of catheter is more difficult and it generates too small lesion size. Moreover, RFA allows the use of expandable electrodes and placement of multiple probes, which is not possible with cryoablation.

Laser Ablation

Percutaneous Laser Ablation (LA) is a technique for image-guided local tumor ablation within solid organs. Laser illumination of the optical fibers introduced into the target lesion delivers energy to the tumor [START_REF] Pompili | [END_REF]], the targeted tissue is heated by the absorbed laser energy and evaporates or sublimates. It can be used on benign and malignant lesions in various organs. The main applications currently involve the destruction of primary and secondary malignant liver lesions. Laser ablation can make deeper lesions and does not predispose to blood clot formation. Yet, there still is a risk of damaging surrounding structures and a very limited experience with laser even during open-heart surgery. 8

Electroporation

Electroporation is a molecular biology technique in which an electrical field is applied to cells to increase the permeability of the cell membrane, that allows chemicals, drugs, or DNA to be introduced into the cell [START_REF] Neumann | [END_REF]]. Irreversible electroporation (IRE) applies short, high-voltage pulses to the tissue to permeabilize cell membranes and induce cellular apoptosis [START_REF] Weiss | [END_REF]]. It can be relevant and efficient in organ like the pancreas, where the other technologies presented previously lack safety and efficacy owing to the risk of thermal damage to adjacent structures and significant thermal dissipation from large adjacent blood vessels. IRE is now used to treat pancreatic cancer previously thought to be unresectable. 9 One potential benefit of IRE is that it generates less heat than RFA. Also, because the heat is generated during the pulses but dissipates while the electrode rests, there is not a constant cooking effect as with RFA, leading to potentially less damage to non targeted tissues. However, clinical experiences are still limited before establishing the efficacy of IRE, whereas RFA have been time-tested for nearly three decades. Multiple needle placements are required within a prescribed distance, which can be challenging, and parallel placement of the probes may be hindered by issues such as intervening ribs. The actual timing of imaging follow-up and the best modality are still being determined. IRE also poses a unique challenge to the current imaging response criteria because unlike other ablation modalities, there is a marked decrease in the size of a successfully treated lesion [Narayanan 2013].

Technical Background

Computational models of RFA are divided into 3 different modeling parts: an electrical heating model, a heat transfer in living tissue model, and a cellular necrosis model. These models are coupled as illustrated on Figure . 2.6 and described in the following sections. [START_REF] Tungjitkusolmun | Three-dimensional finite-element analyses for radio-frequency hepatic tumor ablation[END_REF], Panescu 1995]. The temperature and the voltage are usually set to be constant numbers on the surface of RF electrode probe needle and the potential field generated around the probe is solved numerically, using the Laplace equation (Eq. 2.1) for the electrical field in conductive media [Doss 1982].

Source electrical energy

Heat Propagation in

∇.σ∇V = 0

(2.1) where V is the voltage and σ, the electrical conductivity. The Laplace equation (Eq. 2.1) allows to compute the current density J and the electric field intensity E, using the distributed heat source due to Joule loss q, given by Eq. 2.2:

q = J.E (2.2)
which can also be written as:

q = |∇V | 2 (2.3)
Even if this approach is theoretically accurate, it requires a fine volume mesh on the surface of the very thin probe tips and then suffers from a high computational cost. Furthermore, the exact position of the whole electrode probe is not always available and the imperfect needle positioning has been found to severely affect the outcome of RFA procedure [Khlebnikov 2010]. The electric field strength is not high, most of the Joule effect heating is generated within the space 1 mm around the electrode surface [Panescu 1995] and most of the ablation zone at distances farther from the electrode is created by thermal conduction [Schramm 2006]. Moreover, the particular details of the heating point sources appear to have only a limited effect on the final lesion size [Payne 2011]. Thus, in our framework, the electrical heating is modeled with a Dirichlet boundary conditions fixing the delivered temperature either on a sphere, whose radius is defined pre-operatively by the protocol followed by the clinician if the position of the probe is unknown, or on the probe model itself if we know the probe position.

Biophysical Models of Heat Transfer in Tissue

Computing heat diffusion in biological tissues amounts to solving the coupled bioheat equations derived from the theory of porous media (Eq. 2.4), where each elementary volume is assumed to comprise both tissue and blood with a certain fraction [Nakayama 2008].

(1 -ε)ρc ti ∂T ti ∂t = (1 -ε)Q source +(1 -ε) ∇ • (d∇T ti ) dif f usion +H(T bl -T ti ) (2.4a) ερc bl ( ∂T bl ∂t + v.∇T bl advection ) = εQ + ε∇ • (d∇T bl ) -H(T bl -T ti ) (2.4b)
In theses equations, T , Q, v, ε stand for temperature, source term, blood velocity, blood volume fraction (fraction of blood volume over total volume) and subscripts 'ti' and 'bl' refer to tissue and blood phase respectively. The other parameters are listed in Table 2.1. Different models similar to (Eq. 2.4) have been proposed in the litera- ture [START_REF] Roetzel | [END_REF], Khaled 2003, Nakayama 2001, Shrivastava 2009]. Some models include a metabolic heat generation term [Mahjoob 2009], which can be neglected since it has been shown to be insignificant in thermal ablation [START_REF] Labonté | Numerical model for radio-frequency ablation of the endocardium and its experimental validation[END_REF]]. The advection term in the blood bio-heat equation (Eq. 1.1b) is sometimes simplified by considering an averaged perfused tissue volume [Shrivastava 2009] but we preserve it in our model as it is essential when modeling the heat transfer between the tissue and the small vessels with low blood velocity magnitude. To account for the heat transfer associated with the transcapillary fluid exchange, some models include an additional perfusion heat transfer term [Nakayama 2008]. However, it has been proven that this term can be neglected as the blood in the capillary network usually reaches equilibrium with tissue temperature [START_REF] Chato | Heat transfer to blood vessels[END_REF]]. Some bioheat models examined countercurrent heat transfer in arterial-venous vessels (the venous flow is warmed through heating from the nearby arteries) [Mitchell 1968[START_REF] Roetzel | [END_REF][START_REF] Chato | Heat transfer to blood vessels[END_REF].

One common simplification of the coupled bio-heat equations is the Pennes model [Pennes 1948] where the blood temperature is assumed constant, which is valid within and close to large vessels, where the blood velocity magnitude is high. In this case, the coupled bio-heat equations (Eq. 2.4) reduce to one unique equation with only one temperature distribution T describing the temperature inside the liver:

ρc ti ∂T ∂t = Q + ∇ • (d∇T ) + c bl ω(T b0 -T ) (2.5)
where ω represent the blood volumetric perfusion rate. This simplification has been used widely in the literature to model the electricalthermal heating process happening in RFA. However, this model assumes that the blood is a volumetric heat sink and that it is uniformly distributed throughout the tissue [START_REF] Peng | [END_REF]]. Hence the term accounting for heat-transfer convection between tissue and blood in the Pennes equation is oversimplified. This simplification may be suitable for low vascularized organs like the heart or the kidney and so it can be used more accurately to model cardiac ablation or kidney cryotherapy [Linte 2013b, Talbot 2014] but it is not appropriate for the liver which is highly vascularized. Therefore, we have to account for the perfusion in the parenchyma in another way. Studies [START_REF] Peng | [END_REF]] demonstrate that microvasculature perfusion does not act as a spatially homogeneous heat sink, which invalidates the fundamental assumption of the Pennes model when applied to the liver in which different types of vessels are present. For these reasons another simplification has been proposed, which results in the Wulff-Klinger (WK) model [Klinger 1974, Wulff 1974]. It assumes equilibrium between tissue and blood temperature which is accurate for highly perfused organ with small vessels where the blood velocity magnitude is low and accounts for the directional effect of the blood flow, with a blood volume fraction ε 1:

(

1 -ε)ρc ti ∂T ∂t = (1 -ε)Q + (1 -ε)∇ • (d∇T ) -ερc bl v • ∇T (2.6)
This equation holds in a porous medium where tissue is dominating. The main difference between these two models lies in their cooling terms (last term of the righthand side). The former is a reaction term and acts as a volumetric homogeneous heat sink whilst the latter is an advection term accounting for the directional effect of blood flow on the tissue temperature. It can be shown [START_REF] Peng | [END_REF]] analytically and computationally but in 1D that the coupled bio-heat equations (Eq. 2.4) can be simplified into the Pennes model in the case of large vessels, and to the WK model in the case of small vessels.

Models of Cellular Necrosis

Thermal treatments aim at transporting heat energy within the cancerous tumor then creating a zone of dead cells surrounded by tissue which could eventually recover after the ablation procedure. For a cell to go to an apoptosis state, a critical temperature should be exceeded during a sufficient time duration, due to the thermal tolerance of cells. Various cellular necrosis models [Berjano 2006, Diller 2008, ONeill 2011] have been studied. The simplest one is to use a single temperature threshold above which cells stop instantaneously functioning, and below which cells remain fully functional. Although different values have been used for this threshold, the in vivo lesion volume (i.e. necrotic tissue after RF ablation) can be defined by the volume enclosed by the 50 • C isothermal surface [Panescu 1995], which has been widely used for RFA computational simulations. This isothermal model does not take into account the duration for which cells are at high temperatures although this factor has an influence on the cellular response since tissue damage is a function of both temperature and time [ONeill 2011].

Arrhenius-based models have been proposed [Henriques Jr 1947, Moritz 1947], which extend the Arrhenius law proposed for chemical reaction rates. The rate of cell damage is then proportional to exp(-Ea/RT ), where Ea is an activation energy, R is the universal gas constant, and T is the temperature. In this case, the tissue damage increases linearly with time and hyperbolically with temperature.

Thermal dose models have also been proposed. The earliest work was based on the concept of "equivalent minutes" [Sapareto 1984], where they characterize the thermal dose in an equivalent time at a temperature fixed at 43 • C. The thermal dose in equivalent minutes at 43 • C can be computed at each position according to this empirical model:

t 43 = tc t=0 C (43-T ) ∆t (2.7)
where t 43 is the thermal dose in equivalent minutes at 43 • C, T is the average temperature during the time interval ∆t, and C is a constant resulting from the Arrhenius model:

C = 0.5 when T ≥ 43 • C C = 0.25 when T < 43 • C
A thermal dose of 240 (in equivalent minutes at 43 • C) is considered lethal. This model is based on the concept of a thermal dose, whereby cumulative damage leads to cell death. However, the possibility for a cell to recover is not allowed with this model as the thermal dose is monotonically increasing.

In [Linte 2013b], the authors define cumulative exposure as a measure of induced tissue injury. Cumulative exposure is the area under the temperature-time curve calculated on a voxel-basis over the duration of the ablation, which takes into account the voxel exposure prior to reaching the cell-death temperature (i.e., reversible damage and lesion penumbra) and also the voxel exposure beyond cell-death temperature for 5 s or longer (i.e., irreversible damage and core lesion). Their proposed exposure criterion yields consistent results with Arrhenius and the 50 • C isotherm tissue injury criteria but at a lower computational expense.

Cell-death models based on the Arrhenius law are widely used [START_REF] Pop | A theoretical model for RF ablation of kidney tissue and its experimental validation[END_REF], Chen 2009a, Chang 2004, Schwarzmaier 1998], but are limited since the interpretation of model parameters characterizing the cell damage formulation is not easy to find. These models compute the degree of tissue damage at a tissue location r, using [Pearce 1995]:

Ω(r, t) = ln [C(r, 0)] [C(r, t)] = t 0 Ae E a RT (r, τ ) dτ (2.8)
where Ω(r, t) is the thermal damage index (dimensionless), τ denotes the exposure time (s), A is the frequency constant (s -1 ), E a is the activation energy (J mol -1 ) and R is the universal gas constant (8.314 J mol -1 K -1 ). The thermal damage is defined as the logarithm of the ratio of the original concentration of native tissue at time zero: C(r,0), to the remaining native state tissue at time t: C(r,t). A value of 1 (respectively 4.6) for Ω represents 63% (respectively 99%) of the cells in a damaged state. In some cases, this model is coupled to the heat transfer model: the perfusion term is multiplied by α = 1 Ω , since once tissue coagulation occurs, tissue perfusion ceases [Chang 2004].

Those models are sensitive to small changes in parameters, and they are not able to predict cellular injury over a wide hyperthermic temperature range and throughout the entire heating process [Feng 2011]. To tackle these drawbacks, twocompartment models have been proposed containing either fully alive or fully dead cells [START_REF] Feng | [END_REF], Feng 2011]. In order to characterize the biological state changes, transitions between the two compartments are usually modeled with first-order rate processes.

Finally a three-compartment cell death model [ONeill 2011] can be used to simulate both fast and slow cell death over a temperature range extending to 100 • C. This last model uses a single continuous function where a backward, recovery process is accounted for by adding a vulnerable state of the cell from which cells can either die or heal and return to the alive compartment. In this model, thermal-based necrosis depends on both the heating temperature and the heating time. This model accounts for the recovery process of the cells, allows to define three different states for the cell properties, however, it required to estimate several parameters.

State of the Art of Existing Radiofrequency Ablation Models

Computational Models of RFA

Several research groups are working on methods to extract information from the images, and on models to simulate the RFA behaviour. More specifically, such models could potentially calculate the size, the shape and the location of the necrotic area, given the position and settings of the ablation probe in the abdomen. Thus they could give additional guidance to clinicians in deciding where to place the heating probe and for how long heating must be applied in order to fully ablate the tumor and limit the number of recurrences. For instance, the placement of the probe has been optimized in [START_REF] Schumann | [END_REF]] by using segmentation masks and by the exploration of the set of pareto-efficient solutions, in order to improve the established planning of RFA.

Several approaches have been developed to describe and simulate RFA of liver tumors. They differ in their choice of the biophysical phenomena considered and the type of experimental data used to design and validate them.

All simulations are based on the bio-heat equation considering a cooling effect that is either expanded [Chen 2009b] or localized at neighboring vessels [Altrogge 2012, Jiang 2010, Payne 2011, Audigier 2015a]. Those studies couple radiofrequency electrical fields to thermal transport. They have been developed to compute heat diffusion in the liver, predict the temperature distribution during the procedure, study the biophysical impact of the temperature and eventually evaluate the optimal placement of the RFA probes [Kröger 2010].

Early studies use simplified 3D geometries to determine tissue temperature distribution during RF hepatic ablation [START_REF] Tungjitkusolmun | Three-dimensional finite-element analyses for radio-frequency hepatic tumor ablation[END_REF]]. Other studies [Chen 2009b, Jiang 2010] simulate the heat transfer on generic human anatomies but they do not rely on patient-specific data. Few authors [Altrogge 2012, Audigier 2015a, Payne 2011] have proposed to simulate RFA on realistic subjectspecific geometries extracted from images and only [Audigier 2013a] has personalized biophysical parameters on patient data in order to minimize the discrepancy between simulated and measured necrotic (ablated) regions, which is paramount since parameters in the literature often come from experiments of tissue from different species ex-vivo [Hall 2015].

Some authors simulated the bio-heat equation based on the most common discretization method: the Finite Element Method (FEM) [Chen 2009b, Altrogge 2012, Jiang 2010], Instead we introduced the Lattice-Boltzmann Method (LBM) [Audigier 2015a, Audigier 2013a]. The traditional method to discretize RF ablation PDE is indeed the Finite-Element Method (FEM). However, the use of classic FEM is often more computationally demanding (execution time is usually in the range of hours, which is not suitable for clinical purposes nor therapy guidance) than using recent numerical methods such as the Lattice Boltzmann Method (LBM) [START_REF] Rapaka | LBM-EP: Lattice-Boltzmann Method for Fast Cardiac Electrophysiology Simulation from 3D Images[END_REF]]. Other authors proposed simplified model using a weighted distance-based method [START_REF] Rieder | [END_REF][START_REF] Schumann | [END_REF] to approximate the ablation zone in real-time utilizing the GPU, but the analysis of the effect of using an approximation instead of the complex numerical simulation has not been done yet. Furthermore, a preprocessed thermal equilibrium representation of the liver parenchyma and blood vessels is needed to incorporate its cooling effect. In [START_REF] Rieder | [END_REF]], they only compare the ablation zone given by their method with the coagulative necrosis mask from the numerical simulation in complex vascular situations, indicating a good correlation.

Model Simplification

In most of the developed models, numerous simplifications are assumed, such as the homogeneity and the isotropy of the tissue, a blood perfusion rate unaffected by the heating process, no boiling of tissue during heating [Berjano 2006]. Also the tissue itself is neither isotropic nor homegeneous, particularly in diseased states, which are cases where ablative treatment are commonly used [START_REF] Payne | [END_REF]]. In this work, the properties of the tissue (heat capacity, etc) vary with respect to its state and different parameters are related to different location inside the liver. The reac-tion coefficient, H and the porosity, ε are related to large vessels, and parenchyma respectively. Nevertheless, some models include the water evaporation during the RFA procedure [Yang 2007, Pätz 2010], porous circulation [Crezee 1990], or blood coagulation [Chang 2004]. An important issue is the relationship between tissue characteristics and temperature. RF ablation procedures involves heating at temperature around 105 • C, the parameters values found in the literature may not be correct in this range of temperature. Some models ignore the thermal dependence of the heat diffusivity of tissues [Chen 2009b, Jiang 2010], but it is still not clear if this dependence of the tissue properties have an important impact on the RFA output. For this reason, we keep in this work the definition used in [Payne 2011], where the heat diffusivity is temperature-dependent (Table 2.1).

The liver is a very complex thermal organ due to its inhomegeneity. It is composed of three different types of blood vessels (hepatic arteries, portal veins, hepatic veins), which have different diameters and different flow velocities [START_REF] Tungjitkusolmun | Three-dimensional finite-element analyses for radio-frequency hepatic tumor ablation[END_REF]]. The large blood vessels have a significant cooling effect on RF lesions as it has been shown on ex-vivo porcine liver [Lehmann 2009]. The heterogeneous cooling effect of large neighboring vessels have been studied [Patterson 1998], but the effect of the perfusion is still not totally defined. For these reasons, the vascular system of the liver has to be considered, and we model it in this work but it is often neglected or simplified in most of the studies. Moreover, blood flow circulation is often not computed based on patient-specific clinical information, whereas it is the case in this work. To give an example, the blood temperature in the large vessels is unaffected by the thermal field in the surrounding tissue in [START_REF] Tungjitkusolmun | Three-dimensional finite-element analyses for radio-frequency hepatic tumor ablation[END_REF]]. Other models [Pennes 1948, Panescu 1995[START_REF] Labonté | Numerical model for radio-frequency ablation of the endocardium and its experimental validation[END_REF]] used convective boundary conditions at the tissue -vessels interfaces to simulate the heat transfer due to blood flow, but in order to achieve accurate prediction of RFA output, blood flow needs to be included in the model [Jain 2000].

It is also particularly important to take into account the effect of perfusion, which is often neglected in these models, it is done in this work, as in [Audigier 2015a, Payne 2011] where the cooling effect due to venous flow in the liver parenchyma is considered.

Model Validation

Subject-specific modeling of liver tumor radiofrequency ablation (RFA) needs further validation before impacting clinical use.

Up to now, the comparison between simulated and measured necrotic regions has been used by several authors [Altrogge 2012, Payne 2011, Audigier 2015a] as the main criteria of success in calculating the effect of RFA on abdominal tumors, for either model validation or personalization. In [Payne 2011], the authors compare the simulated volume with the measured one on two pigs and they have an average error of 40% on the volume. Whereas, in our work, point-to-mesh errors of 5 mm on average are computed between the simulated and the actual lesion.

Computational models of RFA depend on a large number of different parameters, which depend mostly on the patient, the temperature or the current state of the tissue. They also have a high computationnal cost. RFA model validation is then challenging with clinical data but is more feasible on pre-clinical data in a more controlled environment. This extensive validation step is important prior to any translation into clinical settings. It is paramount because the size and the location of the RFA lesion has to be predictable and controllable for clinical applications. A comprehensive understanding of detailed multi-physics model is required first to simplify only the aspects that are not needed, which could then enable model personalization or simplification and eventually computational models of RFA can be helpful for interventional guidance or therapy planning. Several experimental validation of RFA modeling have been proposed, but they rely on ex vivo or in vitro data [Linte 2013a, Jain 2000]. In [Linte 2013a], an ex vivo study on two beef muscle samples is performed to validate an endocardial tissue RFA model by comparing the predicted versus the experimentally measured temperature at 2 locations. Validation on in vivo data has been performed as well [Payne 2011], but the validation is done by comparing the volumes of the simulated and post-operative lesions on two pigs only, with nominal biophysical parameters coming from the literature.

Model Personalization

The accuracy of the RFA computation is challenged by the inter-subject variability in anatomies and tissue characteristics [Hall 2015]. Tissue properties depend upon both the tissue state and the temperature. There is little information available on tumor or disease tissue, which are of interest in our case [START_REF] Payne | [END_REF]]. Because tissue properties are patient-specific and can depend on the current state of the tissue, a proper estimation of those parameters is needed but has been often overlooked in the modeling of RFA.

In cardiac modeling, model personalization, i.e. estimation of patient-specific model parameters has been studied. The tissue conductivity as well as the biophysical model parameters for the electrophysiology are estimated in [START_REF] Relan | [END_REF]], whereas in [Marchesseau 2013], the mechanical parameters of an electromechanical model of the heart are personalized.

In RFA modeling, studies [Altrogge 2012, Chen 2009b, Jiang 2010] have investigated the Finite Element Method (FEM) to simulate the heat transfer on generic human anatomies. Simulations with animal-specific [Payne 2011] or patient-specific anatomies were also recently considered [Audigier 2013a] with the inclusion of cooling effects computed from simulated hepatic venous flow and hepatic parenchymal flow. However, nominal tissue parameters were employed in these studies with values often based on ex vivo experiments on animal tissue sometimes with a large varying range between published studies [Hall 2015]. In [Altrogge 2012], the uncertainty associated with biophysical tissue properties (electrical conductivity and thermal diffusivity) is taking into account. The RFA output is evaluated for certain realizations of tissue parameters at an extremely high computational cost. However, they do not tackle the problem of patient-specific parameters by trying to obtain more accurate material parameters, as the uncertainty is considered as an intrinsic attribute of the modeling process.

In this work, we first attempt to estimate patient-specific parameters by minimizing the error between the simulated necrosis and the post-operative registered ground-truth, in the case of multiple lesions ablated per patient [Audigier 2014a]. However, the necrosis of tissue is the resultant of several combined physical phenomena, mainly the heat transfer and cell death mechanisms, meaning that a given ablated region may be explained by several combinations of parameters. In addition to this identifiability issue, the size of the tumor extent can only be known reliably from post-operative imaging and its shape may be highly asymmetric [Hall 2015], which makes it difficult to eventually update the ablation plan during the procedure, which is the clinical end-point. A method that relies also on pre-operative or interventional data for personalization is therefore required for RFA models to be clinically useful. Those observations are complementary information to the necrotic extent for model evaluation and personalization. A patient-specific tool showing the extent of ablation given the probe position, the heat duration, personnalized parameters and patient images will potentially be beneficial in providing a personalized treatment planning and guidance, as it could improve the current clinical outcome.

Surrogate Model for Optimization of the Probe Placement

Inserting the RFA probe in the liver parenchyma requires considerable experience and constitutes a significant mental task for the clinician. During the procedure, multiple criteria for all possible trajectories have to be taken into account like the effect of the large blood vessels, the distance to the lungs, the ribs, etc. For these reasons, several methods to determine the optimal access path have been proposed [Chen 2009a, Altrogge 2012[START_REF] Schumann | [END_REF].

Some automatic and semi-automatic approaches concentrate on calculation and simulation of the ablation zone and the resulting optimal probe placement with respect to the coverage of the tumor. In [Chen 2009a], the probe placement is optimized on a simplified geometry such that the resulting necrotic lesion totally covers the targeted area. The probe placement can be optimized by maximizing the minimum temperature inside the tumor area such that the volume of destroyed tumor tissue is maximized [Altrogge 2012]. Whereas in [START_REF] Schumann | [END_REF]], the optimization of the probe placement is based on numerical optimization and on image processing: segmentations of the area of interest like vessels, lungs, ribs, tumors, liver.

Other methods propose to optimize the planning of the trajectory for percutaneous needle insertions by taking into account the critical structures but also the instrument shape and penetration angle. Methods taking theses constraints into account in order to automatically compute insertion trajectories [START_REF] Villard | [END_REF], Villard 2005, Baegert 2007, Seitel 2011] have been proposed. A simplified radiofrequency ablation simulator can be used [Villard 2003] to get an automatic planning of needle positions considering the tumor environment. The necrosis zone is represented by spheroid deformed by the surrounding vessels to take into account the heat-sink effect. In [Seitel 2011], they use pareto optimality to have a weightindependent approach whereas in [Baegert 2007] the planning of the automatic trajectory in the abdomen is done by a weighting of the different constraints.

MR Thermometry

MRI provides excellent soft-tissue contrast anatomical images as well as a continuous temperature monitoring, thanks to the combination with proton resonance frequency shift (PRFS) based MR thermometry [Ishihara 1995]. Hence, MR thermometry could estimate the resulting lesion sizes of any thermal ablation and thus provide an adequate therapy end point [START_REF] Quesson | [END_REF]]. However, precise temperature measurements using PRF require to elaborate correction techniques, if motion and magnetic field variations are present during the treatment [De Senneville 2007]. More recently, thermal tissue properties have been estimated using MR thermometry during HIFU [Cornelis 2011].

Thermal Tissue Properties Acquisition during HIFU based on MR Thermometry

The development of MR hybrid hyperthermia systems with thermal imaging capability has opened the door to acquisition of thermal tissue properties during heating.

In [Cornelis 2011], the authors evaluate quantitatively in vivo the tissue thermal properties during high-intensity focused ultrasound (HIFU) heating of pig kidneys monitored in real time by volumetric MR thermometry. The analysis of the volumetric temperature distribution during the HIFU intervention allows the determination of the thermal parameters, and may therefore improve the quality of the planning of noninvasive therapy with MR-guided HIFU.

MR Thermometry Filtered by Bioheat Equation

Extending MR temperature imaging to include the entire targeted area and its vicinity with a high accuracy is a difficult problem as the trade-off among spatial resolution, temporal resolution, and volume coverage is inherent in MRI. In addition, the precision of real time MR-thermometry for therapy guidance is generally limited by the available signal-to-noise ratio (SNR) and the influence of physiological noise.

In [De Senneville 2013], they generate 3D temperature maps by combining 2D slices measured by MR thermometry and a priori knowledge of 3-D data derived from predictions based on a physical model (the bio-heat transfer equation model) using an extended Kalman filter (EKF). This is done to monitor HIFU intervention but could be applicable to RFA.

Chapter 3 Based On: the conference paper [Audigier 2013a] and the extended journal version [Audigier 2015a] for the description of the model.

Modeling of Radiofrequency Ablation

Introduction

Turning computational models of RFA into clinical practice remains challenging. The lack of an integrated, efficient, patient-specific framework for RFA modeling based on patient data is the major difficulty encountered. Studies based on in vivo animal images have been reported, but time-consuming FEM computations are still required [Payne 2011]. Because of the expensive computation time needed by FEM, current approaches cannot enable a personalization of the model-based planning of RFA in a clinical setting. Therefore a fast and efficient model would help to optimize the treatment protocol pre-operatively.

As a first step towards efficient patient-specific planning of RFA, we propose an integrated multi-physics approach that combines, the computation of the hepatic blood flow with a biophysical model of the heat transfer and cellular necrosis to simulate RFA therapies in patients and predict the extent of ablation. The computational model of heat transfer and cellular death during RFA is only based on volumetric image of the anatomies estimated from CT (Sec. 3.2.3).

This chapter presents the Lattice Boltzmann Method (LBM) for heat in biological tissue. This method is a kinetic-based approach for fluid flow computations which has been successfully used as an alternative numerical method for solving Navier-Stokes type equations [Chen 1998]. In a first step, the computation of the hepatic blood flow is done using a Computational Fluid Dynamics (CFD) method in the large vessels and Darcy's law in the parenchyma. We will refer to this method throughout the manuscript as the Decoupled FEM Porous -Vessels CFD Model. To go further, an improved hepatic blood flow computation is presented, in which we rely on LBM to compute not only heat diffusion, cellular necrosis as previously proposed but also blood and parenchyma flow in the liver tissue. This latter method is based on a CFD solver which incorporates a porous part to deal with the liver parenchyma. We will refer to this method throughout the manuscript as the Fully Coupled LBM CFD Model. This second framework is particularly efficient for the personalization as it provides a fast solver and naturally accounts for the flow transition between veins and parenchyma. This framework is adapted to situation where no temperature map is available.

A detailed convergence analysis against an analytical solution of the heat transfer equation in Sec. 3.4.1 and a computational efficiency study in Sec. 3.4.2 are presented.

Method

As illustrated in Fig. 3.1, and detailed in the following section, we first estimate a comprehensive level set representation of the liver, including parenchyma, blood vessels and tumors from medical CT images. Then a computational model of heat propagation, cellular necrosis and blood flow through the vessels and liver is solved to estimate the extent of the ablated tissue.

Model of the Patient Hepatic Blood Circulation System

The heat diffusion in liver tissue, involved during the RFA procedure depends on the velocity field of the hepatic blood flow. In this section, we first propose a Decoupled FEM Porous-Vessels CFD Model to estimate the blood flow in the liver parenchyma. The 3D blood flow is first computed on the large vessels and then the porous flow in the liver parenchyma is computed using the Finite Element Method (FEM) on a tetrahedral mesh. As this method relies on FEM, meshing techniques have to be employed, which increases the complexity of the model and the computational errors. For this reason a second Fully Coupled LBM CFD Model is proposed, where the blood flow in the large vessels as well as in the liver parenchyma are computed simultaneously on a Cartesian grid using an Lattice Boltzmann (LBM) implementation.

Model of Heat Diffusion in

Decoupled FEM Porous-Vessels CFD Model

In this section, we present our first approach: the Decoupled FEM Porous-Vessels CFD Model. The modeling pipeline is illustrated on Figure 3.2, and Algo. 1 presents the RFA computation method.

Algorithm 1 Computational model of RFA Update temperature T using Pennes model in the large vessels, Wulff-Klinger model elsewhere

6:

Update cell-state Heat transfer in liver tissue is highly dependent on the blood flow circulation. To solve the WK model (Eq. 2.6), we need the blood velocity field v everywhere in the parenchyma. The blood inside the liver is modeled as a Newtonian fluid with pre-specified density ρ and viscosity µ. v is calculated according to Darcy's

Model of Heat Diffusion in Liver Tissue

Bi-domain bioheat model

CFD in large vessels

Flow in porous media

Cell necrosis

Patient Data

CT Images

Image Processing

Ablation Extent Prediction

Segmentation (Sec 3.2.3) Registration

Probe location

Post-op data Pre-op data

(Sec 3.2.2.1) (Sec 3.2.1) (Sec 3.2.1) (Sec 3.2.2.2) Figure 3.2:
Steps of the proposed method using the Decoupled FEM Porous-Vessels CFD Model (green: input, blue: processes, red: output).

law [Brinkman 1949]:

v = -κ/(µε 2/3 )∇p (3.1)
where p is the pressure. Because of mass conservation, this amounts to solving the Laplace equation:

∇ • (-κ/(µε 2/3 )∇p) = 0 (3.2)
At the border of the liver, Neumann boundary conditions are employed (no flow is leaking). Dirichlet boundary conditions are applied at the tip of the portal and hepatic veins, to define the pressure drop between them. As we cannot estimate these pressures in vivo, we rely on a CFD model of the hepatic venous circulation system to estimate them (Fig. 3.3). We used a full 3D CFD solver (unsteady incompressible Navier-Stokes equations with viscous terms, expressed in an Eulerian framework which embeds the domain boundary using a level set representation of the segmented vessels [Ralovich 2012]). From the segmentation, a tetrahedral multi-domain mesh is generated based on the resulting multi-label mask image using CGAL [Alliez 2011] (www.cgal.org) to compute the porous flow, which is then calculated using FEM on the linear multi-domain tetrahedral mesh. The resulting flow is tri-linearly rasterized on the Cartesian grid after computation. In this work, the effect of heat on the viscosity of the flow is neglected to decouple flow-related from the heat diffusion calculation for computational efficiency. CFD and porous flow are calculated only once at the beginning of the algorithm.

Let ϕ vc in be the vena cava inflow, ϕ p the portal vein inflow and ϕ vc = ϕ vc in + ϕ p the vena cava outflow (conservation of mass, the hepatic artery is neglected in this study but could be added without modification into the framework). We also set the vena cava outlet pressure p 0 = 3 mmHg in the range of physiological values of healthy patients. The values are listed in Table 2.1. First, we compute the 3D blood flow and pressure distribution within the vena cava and hepatic veins as follows: a plug profile velocity field is applied at the inlets (squares in Fig. 3.3), computed from the outflow ϕ p and the cross-sectional area of each inlet. The CFD calculation gives the downstream pressures p - i and the 3D blood flow ϕ i for each inlet of the hepatic vein. Then we estimate the upstream pressure p + , assumed constant, of the portal vein outlets (circles in Fig. 3.3). We solve Darcy's law and optimize over p + such that the computed perfused flow through the hepatic vein inlets ϕ i matches the one computed at the first step using the 3D CFD model. Then once p + is estimated, we compute the blood flow using the 3D CFD solver. Finally, we compute the blood velocity field inside the parenchyma using Eq. 3.1.

Fully Coupled LBM CFD Model

In order to go further and reduce the rasterization errors introduced by the FEM part of the former model, an improved hepatic blood flow computation is presented: a Fully Coupled LBM CFD solver. This method is based on a computational fluid Dynamics (CFD) solver which incorporates a porous part to deal with the liver parenchyma, and is implemented using the Lattice Boltzmann Method (LBM) to compute simultaneously the blood flow in the large vessels and in the liver parenchyma. The modeling pipeline is modified and illustrated in Fig. 3.4.

As detailed in the following section, we still rely on an estimation of a compre-

Model of Heat Diffusion in Liver Tissue

Bi-domain bioheat model CFD in porous media

Cell necrosis

Patient Data

CT Images

Image Processing

Ablation Extent Prediction

Segmentation (Sec 3.2.3) Registration

Probe location

Post-op data Pre-op data

(Sec 3.2.2.1) (Sec 3.2.1) (Sec 3.2.2.2) Figure 3.4:
Steps of the proposed improved method (green: input, blue: processes, red: output). The blood flow is computed in the large vessels and in the parenchyma simultaneously.

hensive level set representation of the liver, including parenchyma, blood vessels and tumor from medical CT images. The computational model of heat propagation and cellular necrosis remains unchanged. However the blood flow computation through the vessels and liver parenchyma is done simultaneously. This method enables to estimate the extent of the ablated tissue as well. Algo. 3 presents the modified RFA computation method, with the differences with the previous model highlighted in red.

Algorithm 3 Computational model of RFA Update temperature T using Pennes model in the large vessels, Wulff-Klinger model elsewhere

5:

Update cell-state

Description

The blood in the main vessels and in the parenchyma are combined in the generalized 3D incompressible Navier-Stokes equation for fluid flow in porous media, thus considerably easing the definition of boundary conditions between vessels and parenchyma, improving on the previous method [Audigier 2013a, Audigier 2015a]. More precisely, writing v as the blood velocity and p the pressure inside the liver, we solve:

∂v ∂t + v.∇v = - 1 ρ ∇p + µ ρ ∇ 2 .v + F F = -µ(1-ε) 2 ρα 2 ε 2 v (3.3)
The added force F represents the total body force due to the presence of a porous medium [START_REF] Guo | [END_REF]]. F depends on the porosity coefficient ε (fraction of blood volume over the total volume) whose default values are 1 in the CT-visible vessels, 0.1 in the porous parenchyma [Payne 2011], and 0.04 in the vessel walls to model an impermeable medium. Experiments have been performed to obtain a sufficiently small porosity (0.04) to avoid flowing through the vessel wall. In Eq. 3.3, ρ stands for the blood density, µ, the dynamic viscosity of the blood and α 2 an effective parameter. At the border of the liver, no flux boundary conditions are used whereas Dirichlet boundary conditions are applied at the inlets of portal vein and vena cava and at the outlet of the vena cava: the portal vein and vena cava inflow, ϕ p and ϕ vc in are fixed as well as the vena cava outlet pressure p 0 (see Figure 3.5 for more details). This method makes the boundary conditions simple to treat: no boundary conditions are fixed on the extremities of the vessels inside the parenchyma thanks to the use of the porosity map. This framework mainly avoids the occurrence of shear stress on the vessel walls due to their much lower value of porosity. Link between the Decoupled FEM Porous-Vessels CFD and the Fully Coupled LBM CFD Models In the Decoupled FEM Porous -Vessels CFD Model(Section 3.2.1.1), the blood flow is computed according to Darcy's law: Eq. 3.1 [Brinkman 1949]. In the Fully Coupled LBM CFD Model (Section 3.2.1.2), the blood flow is computed according to the generalized 3D incompressible Navier-Stokes equation for fluid flow in porous media: Eq. 3.3. At steady state, Eq. 3.3 gives the blood velocity as a function of the gradient of the pressure:

v = - α 2 ε 2 µ(1 -ε) 2 ∇p (3.4) α 2 = d 2 p /150
, where d p is the diameter of the solid particle [START_REF] Guo | [END_REF]]. From (Eq. 3.1) and (Eq. 3.4), the same law should link the velocity and the gradient of the pressure. To this end, the value of d p is set to d p = 1.5mm after computing both coefficients with biophysical values of parameters from Table . 2.1.

Numerical Resolution using LBM

Eq. 3.3 is solved using the Lattice Boltzmann Method (LBM) for fast computation on general purpose graphics processing units (GPU). LBM has been developed for CFD and is now a well-established discretization method [Chen 1998]. In RFA, it has been validated in [Audigier 2013a] through a comparison with an analytical solution, for a similar accuracy as FEM. In this case, LBM is used to compute the porous and blood flow circulation in the liver. To this end, an isotropic Cartesian grid with 19-connectivity topology (D3Q19 scheme defined in Figure 3 A Multiple-Relaxation-Time (MRT) model is used for increased stability [START_REF] Pan | An evaluation of lattice Boltzmann schemes for porous medium flow simulation[END_REF]]. At position x for the edge e i , the governing equation is:

f (x + e i ∆x, t + ∆t) = f (x, t) + A[f eq (x, t) -f (x, t)] + ∆tg(x, t).
(3.5)

In this equation, f (x) = {f i (x)} i=1..19 is the vector of distribution function with f i (x) being the probability of finding a particle travelling along the edge e i of the node x at a given time; c = ∆x/∆t; c 2 s = 1/3; ∆x is the spacing;

f eq i (x, t) = ω i ρ[1 + e i .v cc 2 s ] g i (x, t) = ω i ρ e i .F c 2 s (3.6) ω = {ω i } i=1.
.19 is the vector of weighting factors and A the MRT matrix. The fluid mass density and velocity are computed from the LBM distributions as:

ρ = 19 i=1 f i (x, t) ρv = 19 i=1 e i f i (x, t) + ∆t 2 ρF (3.7)
The fluid mass density and the velocity are updated at every node of the grid for every timestep ∆t.

Differences between the Decoupled FEM Porous-Vessels CFD Model and the Fully Coupled LBM CFD Model

We first proposed the Decoupled FEM Porous-Vessels CFD Model to compute the blood flow in the hepatic vessels and in the liver parenchyma. This model requires two steps: first the computation of the blood flow in the large vessels using a Computational Fluid Dynamics (CFD) solver and then the computation of the blood flow in the parenchyma using Laplace equation and Darcy's law. This method requires mesh generation and a significant computational time. For these reasons, we proposed a second method where both flows (in large vessels and parenchyma) are computed simultaneously on the Cartesian grid using the Lattice Boltzmann Method (LBM).

There are differences in practice between those two methods. The Decoupled FEM Porous-Vessels CFD Model requires two decoupled computations. Those two computations are computationally demanding and we have to fit the values found by the second computation at the tips of the hepatic and portal veins to those found by the first computation. Moreover, the computation of the porous flow is done using the Finite Element Method (FEM). The computed flow is tri-linearly rasterized on the Cartesian grid to be used in the heat transfer model, which introduced interpolation errors. On the other hand, the Fully Coupled LBM CFD Model considerably eases the definition of boundary conditions between vessels and parenchyma: no boundary conditions are fixed on the extremities of the vessels inside the parenchyma thanks to the use of a porosity map, which is generated from the vessel mask. Secondly, no advanced meshing techniques are required, reducing the rasterization errors. In both cases, only the flow going through the portal inlet circled in white are shown. We can see differences between the results given by both computations. The area circled in black shows flow lines, which appear with the first computation and not the second one.

Fully Coupled Model of Heat Transfer in Liver Tissue and Cellular Necrosis

Model of Heat Transfer in Liver Tissue

This model describes how the heat flows from the probe through the liver while accounting for the cooling effect of the main vessels and parenchyma.

As current imaging techniques do not support an accurate measurement of the ratio between blood and liver tissue (large vessels can be clearly identified, but small capillaries are difficult to image), we use two simplifications of the coupled bioheat equation (Eq. 2.4) in the parenchyma and in visible blood vessels. Both equations can be easily implemented in a modular way to cope with tissue inhomogeneity. Continuity between the two models has to be ensured at the extremities of the large vessels. To that end, reaction coefficient (H or R) are chosen large enough in order to reduce the effect of the high temperatures relative to normal body temperature, i.e 37 • C. In our framework, a two-compartment model is used according to the spatial location within the anatomy: either the Pennes model or the WK model is used. Assuming that large blood vessels and surrounding tissue are isolated from each other, we compute the temperature by solving the diffusion equation:

ρc ti ∂T ∂t = Q + ∇ • (d∇T ) (3.8)
everywhere in the domain, to which we add the cooling term: H(T b0 -T )/(1 -ε) when a point belongs to a large vessel, where blood velocity is high, (Pennes model) or -ερc bl v • ∇T /(1 -ε) when it belongs to the parenchyma, where tissue is dominating (WK model). T , Q, v and T b0 stand for temperature, source term, blood velocity and the mean temperature (assumed constant) of the blood in large vessels.

In the case of a Fully Coupled LBM CFD Model, as the porosity ε vary spatially inside the liver, the temperature T is computed by solving the following equations: a reaction-diffusion equation in the large vessels, and an advection-diffusion equation in the liver parenchyma:

ρc ti ∂T ∂t = Q + ∇ • (d∇T ) +R(T b0 -T ) -α v ρ b c b v • ∇T (3.9)
where ρ, c, d are the density, heat capacity, conductivity; subscripts ti and b stand for tissue and blood respectively. Q is the source term, α v , the advection coefficient, R, the reaction coefficient and T b0 the blood temperature (assumed constant) in large vessels. In order to be consistent with the previous model R is set to

R = 27.1 × 10 4 W (m 3 K) -1 and α v = 0.11.
It can be shown [START_REF] Peng | [END_REF]] analytically and computationally but in 1D that the coupled bio-heat equations (Eq. 2.4) can be simplified into the Pennes model in the case of large vessels, and to the WK model in the case of small vessels.

Our model includes the heat sink of all hepatic vessels (veins and arteries) as well as the effect of the blood flow within the parenchyma considered as a porous medium. A weakly coupling model is considered: the blood flow has an influence on the temperature distribution through the WK model and through the reaction term in the Pennes model but the temperature does not affect the blood flow (coagulation is not considered here), which allows us to speed up the calculations since the blood flow distribution is computed only once, at the beginning of the simulation. The CFD solver is run until a steady state is reached.

Cellular Necrosis Model

General Formulation

Tissue necrosis is calculated based on the computed temperatures using a threestate model [ONeill 2011, Breen 2002], coupled with the bio-heat equation. Each cell has a probability to be either undamaged (U), vulnerable (V) or necrotic (N). Those probabilities vary with the temperature spatially and temporally according to the following transition diagram (Eq. 3.10):

[U] β(T ) ---→ ← - γ [V] δ(T ) ---→ [N] (3.10)
where β(T ) = βe T /T k , δ(T ) = αβ(T ) are the rates of cell damage and γ is the rate of recovery. Unlike in [ONeill 2011], a constant α is introduced. This diagram results in three coupled ODEs (Eq. 3.11), solved with a first order explicit scheme on the same grid and with the same time step as the bio-heat equation:

             dU dt = -β(T )U + γV dV dt = β(T )U -(γ + δ(T ))V dN dt = δ(T )V (3.11)
The following property is also imposed at each point:

U + V + N = 1 (3.12)
The initial conditions are chosen as in [ONeill 2011]: U = 0.99, V = 0.01 and N = 0. We solve the resulting three coupled ODEs with a first order explicit scheme at each vertex of a Cartesian grid, yielding a spatially-varying cell state field used in the bio-heat solver. The computation is done on the same grid and with the same time step as the bio-heat equation. For both heat transfer and cell death models, parameters are initially set to values from the literature [Pennes 1948] reported in Table 2.1. The cell death model is strongly coupled to the bio-heat equation as the heat capacity depends on the state of the cell (c U t , c V t and c N t correspond respectively to the heat capacity of undamaged, vulnerable and necrotic tissue) and the conductivity d t depends on the temperature through d t = dt * (1 + 1.61 * (T -310).10 -3 ) as in [Payne 2011]. All the remaining parameters are constant.

Nominal Formulation from [ONeill 2011]

To get the same equation as in [ONeill 2011], we set α = 1, and the two damage rate coefficient are equal. In this case, we rename the parameters for the sake of clarity and Eq. 3.10 becomes Eq. 3.13:

U k f (T ) ---→ ← - k b V k f (T ) ---→ N (3.13) k f (T ) = kf e T /T k (1-U )
and k b are the rates of cell damage and recovery respectively. kf is a scaling constant and T k is a parameter that sets the rate of the exponential increase with respect to the temperature. This diagram results in three coupled ODEs (Eq. 3.14):

             dU dt = -k f (T )U + k b V dV dt = k f (T )U -(k b + k f (T ))V dN dt = k f (T )N (3.14)
In this case, a single damage process incorporates all physiological damage mechanisms, thus the transitions U -→ V and V -→ N are the same. The vulnerable state is an arbitrary position representing the "point of no return", it is not a change in Figure 3.8 represents the solution of (Eq. 3.14) at one vertex over time if a constant temperature of 105 • C is applied.

Numerical Resolution using LBM

We rely on the Lattice Boltzmann Method (LBM) where a statistical description of the system is used to compute heat diffusion and cellular necrosis in the liver tissue. LBM is a new computational method, which discretizes the velocity space on a grid, with mass, momentum and energy conservation conditions. To meet these three conservation conditions, two separate distribution functions are usually used [START_REF] Mai | A thermal lattice Boltzmann model for flows with viscous heat dissipation[END_REF]]. In the Decoupled FEM Porous-Vessels CFD approach, FEM solvers (CFD solver in the large vessels and porous solver in the parenchyma) give the blood flow distribution in the liver and enforce the mass and momentum conservation. Energy conservation is modeled using a distribution function for the thermal energy. Whereas, in the case of Fully Coupled LBM CFD solver, to meet the three conservation conditions of mass, momentum and energy, two separate distribution functions are used as in [START_REF] Mai | A thermal lattice Boltzmann model for flows with viscous heat dissipation[END_REF]]. In this approach, two LBM solvers (one CFD solver in the large vessels and in the parenchyma and one heat transfer solver) are used, the first one gives the blood flow distribution in the liver and the second one the temperature and cell state distributions.

LBM has emerged as a powerful technique for efficient computation of second order elliptic partial-differential equations [Yoshida 2010]. It is used in Computational Fluid Dynamics to solve Navier-Stokes equations [Chen 1998], in heat transfer problem [Dawson 1993], or in cardiac electrophysiology [START_REF] Rapaka | LBM-EP: Lattice-Boltzmann Method for Fast Cardiac Electrophysiology Simulation from 3D Images[END_REF]] to solve reaction-diffusion equations. Contrary to FEM, this discretization method uses a grid, allowing an easy parallelisation. We used an isotropic Cartesian grid and, LBM is performed with a D3Q7 scheme (DnQm denotes m discrete velocities in n dimensions): 6 directions are considered + the current point [START_REF] Rapaka | LBM-EP: Lattice-Boltzmann Method for Fast Cardiac Electrophysiology Simulation from 3D Images[END_REF]] as described in Figure 3.9.

Figure 3.9: The D3Q7 scheme used in the heat diffusion solver from [Yoshida 2010].

Numerically, LBM uses fictitious particles, which perform consecutive collision and streaming processes over the discrete lattice mesh. When there is a non-zero probability that particle distributions move to the same node from different directions, then the probability of having particle distributions at that node with a given velocity direction is changed due to the application of a collision operator. The governing equation at position x = (x, y, z) for the direction defined by the vector e i is given by the two equations (Eq. 3.15, Eq. 3.16).

f (x) = {f i (x)} i=1..7 is the vector of the temperature energy distribution function with f i (x), the probability of finding a particle traveling along the vector e i of the node x at a given time.

f (x) describes the post-collision distribution.

ω = {ω i } i=1..7 = (1/4, 1/8, 1/8, 1/8, 1/8, 1/8, 1/8) is the vector of weighting factors [START_REF] Rapaka | LBM-EP: Lattice-Boltzmann Method for Fast Cardiac Electrophysiology Simulation from 3D Images[END_REF]] that depends on the lattice connectivity.

First, the collision step is:

f (x, t + ∆t) = f (x, t) + M -1 Ŝ[Mf eq (x, t) -Mf (x, t)] + ω∆tH(T b0 -T (x, t)) (3.15)
and then, the streaming step is:

f i (x + e i ∆x, t + ∆t) = fi (x, t + ∆t) ∀i ∈ 1...7
(3.16) with:

f eq i (x, t) = ω i T (x, t)[1 + e i .v cc 2 s ]
We denote c = ∆x/∆t, c 2 s = 1/4 and ∆x as the spacing. The set of vectors e i is defined as:

[e 1 , e 2 , e 3 , e 4 , e 5 , e 6 , e 7 ] T =   0 1 -1 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 1 -1   At each time step, the entire 3D lattice domain is traversed and for each cell, new distribution function values are computed from its six neighbors (plus itself). The temperature corresponding to the 0 th moment is computed as:

T (x, t) = 7 i=1 f i (x, t) (3.17)
Only the data of the cell and its six neighbors are used, so the lattice can be traversed in any order since values from the neighbors are computed from the previous time step. The LBM offers high parallel scalability, second order accuracy in space and the simplicity of an implementation on a uniform Cartesian grid [START_REF] Rapaka | LBM-EP: Lattice-Boltzmann Method for Fast Cardiac Electrophysiology Simulation from 3D Images[END_REF]]. The relation between LBM equations (Eq. 3.15, Eq. 3.16) and the continuous model can be derived by employing the Chapman-Enskog expansion, (a formal multi-scale expansion) [Chen 1998]. The derivation are given in Appendix A. For stability reason, we use a Multiple-Relaxation-Time model since the numerical stability can be improved by separating the relaxation rates of the conserved and non-conserved moments [START_REF] Humières | Multiple-relaxation-time lattice Boltzmann models in three dimensions[END_REF]. First, f is brought to a new basis in which each component corresponds to a certain moment of the vector Mf (0 th order is the temperature T ). Then each component relaxes to the equilibrium Mf eq with a different relaxation coefficient. Finally, the vector is projected back onto the original seven-dimension space [Yoshida 2010]. In (Eq. 3.15), instead of writing A = M -1 ŜM = 1/τ I, the relaxation towards equilibrium is performed in the moment space, where

M =            1 1 1 1 1 1 1 0 1 -1 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 1 -1 6 -1 -1 -1 -1 -1 -1 0 2 2 -1 -1 -1 -1 0 0 0 1 1 -1 -1            and Ŝ = diag(1, 1/τ, 1/τ, 1/τ, 1/τ 1 , 1/τ 1 , 1/τ 1 )
is the collision matrix in moment space. The relaxation time τ is directly related to the diffusion coefficient D through τ = 1/2 + 4D∆t/∆x 2 [Yoshida 2010]. After a stability analysis, we chose τ 1 = 1.33 [START_REF] Rapaka | LBM-EP: Lattice-Boltzmann Method for Fast Cardiac Electrophysiology Simulation from 3D Images[END_REF]] to get a stable and well-behaved solution. Based on the Chapman-Enskog expansion, this value does not change the accuracy of the solution but it affects its stability.

No-flux boundary conditions are applied at the border of the liver. The boundaries are treated according to the level set representation of the liver computed from the segmentation using linear interpolation without requiring advanced meshing techniques. We use a second-order accuracy model for curved walls presented by [Bouzidi 2001] who proposed a simple boundary condition based on interpolation and the bounce back scheme. This method needs to treat the boundary conditions separately for two cases: ∆ ≤ 1 2 or ∆ > 1 2 , i.e either the boundary is closer to the lattice point which is inside the domain or it is closer to the lattice point which is not in the domain. ∆ represents the fraction of an intersected link in the domain of interest (Figure 3.10 illustrated the description of ∆) and is computed based on the level-set representation of the liver.

Δx

ΔΔx

Boundary wall out 0 -level set If the node from which the post collision values travels (x -e α ∆x) is outside the domain, then the value of f α (x, t + ∆t) is:

for ∆ > 1 2 :f α (x, t + ∆t) = 1 2∆ fᾱ (x, t + ∆t) + 2∆ -1 2∆ fα (x, t + ∆t) for ∆ ≤ 1 2 :f α (x, t + ∆t) = 2∆ fᾱ (x, t + ∆t) + (1 -2∆) fᾱ (x -e ᾱ∆x, t + ∆t)
where e ᾱ = -e α . Finally, we model the heat source term through a Dirichlet boundary condition at the location of the probe. We emulated the RFA protocol by computing a sphere centered at the center of the tumor and with the radius defined pre-operatively by the clinician given the size and the location of the tumor, as illustrated on Figure 3.11.

The temperature of the points inside the sphere is enforced at 105 • C for a duration defined also pre-operatively.

Sphere simulating the probe Tumor RFA Probe

Estimation of Anatomical Model from Patient Images

The pre-and post-operative images are semi-automatically segmented, yielding a detailed anatomical model of the patients liver. The user provides different strokes for the background and the foreground of the image, and the strokes can be modified interactively [Criminisi 2008]. Finally the resulting segmentation is manually refined using itk-snap [Yushkevich 2006]. We generate volumetric binary images of the parenchyma, tumors, hepatic veins, vena cava, portal vein, hepatic artery whenever visible (Figure 3.12). The duration of the full segmentation process depends on the extent of the visible vessels. For example, typically for one patient, it took 20 min to segment the liver mask, 10 min for the tumor, 35 min for the vena cava and 15 min for the portal vein, i.e 1 hour 20 minutes in total. A multi-label mask image is created to identify the different structures of this detailed anatomical model of patient's liver and circulation tree as well. To define the computational domain, a level set representation of the liver, without tumor and vessels is computed.

The Decoupled FEM-LBM CFD solver rely on FEM. In this case, a smooth polygonal surface mesh is created for each region. A multi-label mask image is also created to identify the structures of interest for the computation.

Contrarily, the Fully Coupled LBM CFD solver does not rely on polygonal meshes at all, but the porosity needs to be defined everywhere in the computational domain. In this case, we do not generate a surface mesh but instead a porosity map is necessary. From the vessel masks, a porosity map is created to identify the porous parenchyma and the vessels. Because non-visible, the walls of the vessels are extrapolated using 26 connectivities dilatation of the vessels masks. The vessel extremities, which do not have walls are manually identified.

Heating Power and Cooling Temperature Computation

During the heating period, a Dirichlet boundary conditions is used to fix the delivered temperature either on a sphere, whose radius is defined pre-operatively by the protocol followed by the clinician if the position of the probe is unknown, or on the probe model itself if we know the probe position. In this case, RFA is simulated by imposing as input the temperatures measured at the five thermistors and the four remaining tips temperature are linearly interpolated from these values. The heating stops at time t = t a . During this period (t < t a ) the heating power can be computed, whereas the cooling temperature can be computed when (t > t a ).

Heating Stage

We assume that the measured power is strongly correlated (proportional) to the heat power P (t) delivered through radio-frequency to heat the liver tissue. Proportionality is assumed to account for power dissipation due to electrical resistance, and the unknown surface ratio of the probe being heated. The heat power P (t) delivered to the tissue can be computed at each time step of the simulation from the bio-heat equation according to Fourier's law:

P (t) = α S d t ∂T (t) ∂n dS (3.18)
where S is the probe surface, n is the outer normal at that surface and α, the proportionality coefficient (α = 0.4469), found by matching the peak value of the measured power for pig 1 with the peak value of the simulated power with personalized parameters. This value is then used in all the computations.

Cooling Stage

In the absence of any delivered heat power, the nine tips of the probe cool down at a speed which depends on the local conductivity d t and the heat capacity c t . Thus, during the cooling period t c -t a (cooling stops at time t = t c ), the cooling temperature can be simulated.

Model Verification

The computer hardware used in all the experiments was a Windows 7 desktop machine (Intel Xeon, 2.80 GHz, 45GB RAM, 24 CPUs) with a Nvidia Quadro 6000 1.7 GB with 448 CUDA cores. See text for details.

Experimental protocol

To evaluate the LBM heat transfer solver, its behavior is compared on a regular cuboid domain with an analytical solution. For a source released at x 0 at time t 0 , the 3D analytical solution of the advection-diffusion equation:

∂T ∂t +v • ∇T = ∇ • (D∇T ) is: T (x, t) = M [4π(t -t 0 )D] 3/2 exp - x -x 0 -(t -t 0 )v 2 4D(t -t 0 )
To have the same conditions as the RFA computation on a patient's liver, parameters were chosen to get the heat diffusion in a physiological range. We set a Gaussianshape source of 70 • C at the center of the cuboid at time t = 0 (Figure 3.13) with this set of parameters: D = 0.1 mm 2 /s, v = (2, 0, 0) mm/s, M = 35000 • C.mm 3 , t 0 = -50 s. We initialized the temperature values at each point of the domain with the analytical solution at time t = 0, then the advection-diffusion equation is solved using our LBM solver. The temperature at a typical point of the domain is reported and compared with the analytical solution values. Neumann boundary conditions were used at the border of the domain, which was chosen to be large enough to get rid of the boundary effect at the probed points.

LBM Convergence Analysis

We performed a spatial and a temporal convergence analysis of the solution computed by the LBM solver. The solutions are compared to the analytical solution at one typical point of the domain (Figure 3.14). This analysis confirms qualitatively and quantitatively the accuracy of the implementation of the heat transfer model with LBM.

Temporal Convergence Analysis For a given resolution, an upper and lower bound for the time step were provided by the simulated physics from the Table 3.1: The Root-Mean-Square (RMS) errors between the computed solution and the analytical with respect to the resolution for a fixed time step of 0.01 s .

Resolution RMS 5 mm 11.86 • C 2 mm 2.42 • C 1 mm 0.55 • C 0.5 mm 0.11 • C
LBM implementation and the Courant-Friedrichs-Lewy conditions: |v∆t/∆x| ≤ 1 [START_REF] Courant | Über die partiellen Differenzengleichungen der mathematischen Physik[END_REF]]. A time-step of 75 ms and a resolution of 1 -2 mm appeared to be a good compromise between accuracy and computational cost (Figure 3.14, Right).

Computational Efficiency

A GPU-based version of our model has been developed, relying on the CUDA toolkit (dedicated software for NVIDIA's GPUs) since LBM is easily parallelizable. In our model, the values of interest are the temperature related to the distribution function and the state (alive, vulnerable, or dead) of the cell computed at each vertex separately. This implementation uses classic parallelization methods, similarly to a CPU parallel approach: each thread is dedicated to one vertex and computes the contribution of the temperature and state of the cell for this vertex ensuring a tiled access in memory. Two distribution functions are actually needed (at time t and time t -1: f past and f present ), to avoid concurrency reading and writing when one thread writes a vertex value while another thread attempts to read it. The implementation is described in Algorithm 4.

Algorithm 4 Implementation of LBM RFA Enforce temperature at the probe points.

5:

Compute the corresponding distribution functions f . 6: Compute f and then f 7: Compute the cell state and the temperature T. 8: Update τ , c ti and d. 9: Do f past = f present .

In order to show the benefit of our GPU-based approach, it is compared against a CPU implementation with multithreading. The same synthetic set-up as described in Figure 3.13 is used. As reported in Figure 3.15, experiments showed a maximum speed-up of 11 with multithreading and 22 threads (OpenMP) and 45 with graphical processing units (GPU) implemented in CUDA with respect to a singlecore implementation of LBM. Moreover, after a quantitative verification of the FEM computation against the analytical solution, experiments showed that a 60× speedup was obtained with respect to an FEM implementation on CPU for a similar accuracy.

Computation

Discussion

Summary

We have presented two multi-physics models for efficient patient-specific planning of RFA based on medical CT images. In a first approach, we rely on LBM to solve the bioheat equations and the level set representations of the structures are directly computed from images, but this framework still requires advanced meshing techniques to compute the flow in porous media, adding one step in the processing pipeline with the associated computation time, errors and parameters. To solve this issue, we propose a second model, where we rely on LBM not only to solve the bioheat equations and the cellular necrosis model, but also to compute the blood flow simultaneously in the large visible vessels and in the parenchyma. The coupled computation of the porous and hepatic flow eliminates the difficulties in the setting of boundary conditions in the parenchyma, and the occurence of shear stress on the vessels wall is avoided as the Pennes Model is used in the big vessels where the flow is not accounted for.

Model Limitations

In this study, a two-compartment model was used to describe the heat propagation, coming from two simplifications of the coupled bio-heat equations which are accurate in the cases of small and large vessels respectively, whereas previous studies suggested that heat dissipation may arise precisely in the medium vessels [El-Serag 2003]. This simplification was motivated by the fact that the liver is highly vascularized with different types of vessels of different diameters. For example, the large hepatic veins have a diameter of 6 mm on average and the capillaries, a diameter of 8 10 -3 mm [START_REF] Peng | [END_REF]]. With the resolution of the images used in clinical settings, we cannot segment such small vessels and thus we cannot include them in the current modeling pipeline. Modeling the impact of all vessels is not practical in clinical practice. However, it will be interesting to study the effect of the vasculature on the necrosis, using a synthetic case with a deeply detailed vessel tree for example.

Those experiments are necessary to understand the impact of every hepatic blood vessels. It will be interesting to know until which resolution it is necessary to go in order to get an accurate prediction of the ablation extent. Only then computational models could be simplified and used in clinical settings.

Conclusion

Turning computational models of RFA into clinical practice is a challenging task, which can be beneficial for therapy understanding. Through the use of the Lattice Boltzmann Method, our radiofrequency ablation (RFA) model allows near real time and state of the art computations of heat transfer that are suitable for model-based therapy planning or guidance in the future, even if the target is to go beyond realtime, as we need to stay under 1-2 minutes of computation time for clinical use. Based on our approach, the clinician could test different therapeutical strategies in-silico, assess their outcome before the intervention and finally choose the most appropriate therapy for a specific patient. The model has been verified against an analytical solution, evaluation on clinical data is presented in the next chapter.

Introduction

In this chapter, we describe the application of the framework described in the previous chapter, on fourteen tumors from ten patients. The Decoupled FEM Porous-Vessels CFD Model is used to compute the blood flow (the work presented here has been realized at the beginning of this work and at that time, the Fully Coupled LBM CFD Model was not introduced yet), the nominal formulation from [ONeill 2011] of the cell death model (with α = 1) is used, and different hypotheses are tested in Section 4.4. The comparison with the real outcome extracted from post-operative data in terms of necrotic area suggests a promising correlation between the predicted and actual ablation extent, but also the importance of considering the blood perfusion inside the parenchyma. The generalization of the biophysical model with respect to tissue properties (i.e. we simulated a RFA procedure by using patientspecific geometry and boundary conditions but generic tissue parameters) shows that reasonable results can be achieved (average values of point-to-mesh distance: 10.17±8.52 mm, DICE: 41.8%, sensitivity: 66.9%, positive predictive value: 38.3%). This study presents the generalization of the biophysical model with respect to tissue properties tested on ten patients data, and an analysis of the perfusion effect on the necrotic area. A first experiment of tissue parameters fitting has been performed on one patient. It highlighted a more accurate prediction power but also the need for personalization.

Available Clinical Data

Our database is made of CT scans from ten patients who underwent RFA of one or several hepatic tumors. Figure 6.1 illustrates the case of patient 04. For each patient, pre-and post-operative late venous phase CT are available. We generate volumetric binary images of the parenchyma, tumors, hepatic veins, vena cava, portal vein, without the hepatic artery since only single-phase CT images are available. The clinical data were processed off-line. The results were not available to the interventionist during the RFA procedure. This set of patient data come from the Johns Hopkins Hospital in Baltimore, MD, USA. The computation using the RFA model has been described in the previous chapter and is now evaluated on clinical, retrospective data from ten patients, with fourteen ablations (some patient had several tumors ablated) for whom pre-and post-operative CT images were available. In this chapter, the nominal formulation from [ONeill 2011] of the cell death model (with α = 1) is used to compute necrotic extent and the Decoupled FEM Porous-Vessels CFD Model is used to compute the blood flow. For all patients, nominal tissue parameters were employed (Table 2.1).

Tumor Post-op lesion

Clinical RFA protocol requires that the probe is deployed within the tumor with a probe diameter defined pre-operatively according to the size of the tumor, and then maintained for 7 minutes after the target temperature of 105 • C was reached, as measured by the probe thermistors. For large tumors, the process was iterated with sequentially increasing diameters. In all cases, a single probe placement was utilized (no separate overlapping ablations). After anatomical model extraction, we emulated the RFA protocol by placing the virtual probe at the center of the tumor. Cells around the probe tip within the probe diameter sphere were heated at 105 • C during 7 minutes or twice 7 minutes. In all cases, the computation continued for 3 additional minutes without the probe so that each cell reaches a steady state. In order to evaluate the results of the computation, we compared the computed necrotic region with the patient-specific ground truth. For each patient, the lesion is manually segmented by an expert on the post-operative image and then non-rigidly registered to the pre-operative image. The elastic registration from the Advanced Normalization tools (ANTS) is used [Avants 2009] with the vessels and parenchyma binary images employed as landmarks. Indeed, the thermal induced lesion on the post-operative image creates registration issues, especially in the area of interest (Figure 4.2). Finally the accuracy of the registration is visually checked by the expert.

Computation Time

One minute of ablation is computed in almost one minute. In comparison, using FEM on CPUs, the identical process takes around one hour. To the best of our knowledge, this is the first time that near real-time physiological computations of RF ablation has been achieved. A single probe ablation of 7 minutes is computed in around 8 minutes depending on the liver size, while an increased probe diameter ablation of 14 minutes is computed in around 15 minutes.

Preoperative image

Registered postoperative image Overlay between preoperative and registered postoperative images The shape of the ablated area is also dependent on the heat advection through the liver parenchyma ( Given V m and V s , the volume of measured (respectively simulated) necrotic area, then the DICE score (DICE), the sensitivity (S) and the Positive Predictive Value (P P V ) are defined as: More importantly, in all but one case (patient 01) the computation predicted that the selected protocols completely covered the entire tumor, which is the clinical criterion for ablation planning. The sensitivity is also reasonable (67% on average). The average Positive Predictive Value (PPV) of 38% is low as the values are drastically different for each patient (from 95.6% for patient 02 to 11.9% for patient 09). Some cases presented a computed lesion far from the registered post-operative one. For example, as one can see on Figure 4.12, the computed lesion in red is different and far from the registered post-operative lesion in white for patient 03 (PPV of 14.5%). This may be due to the uncertainty of the actual probe position used in the clinical protocol or to registration errors as discussed in the following section (Section 4.5). Some cases presented a larger necrosis area compared to the ground truth (Figure 4.3) or a smaller one (Figure 4.14,Left). The diffusion coefficient used from the literature was either too high or too low to get a perfect match, as exposed in the following experiment. It may be due to the general state of the liver, which can be cirrhotic or hyperperfused [Van Beers 2001].

DICE = 2|V m ∩ V s | |V m | + |V s | S = |V m ∩ V s | |V m | P P V = |V m ∩ V s | |V s |

Effect of Advection

In order to better understand the effect of advection, we perform the computations on the same data and we remove the blood flow in the parenchyma. The evaluation results are shown on Figure 4.6 and the values are reported in Table 4.2. Qualitatively, we can observe that the advection has an impact on the shape, the extent and the size of the lesion (Figure 4.3 and Figure 4.5). As the blood flow in the parenchyma tends to go from the extremities of the portal vein to the hepatic veins (the extremities of the vena cava), the temperature follows the same path. Notably, if the tumor is close to the portal vein, the advection will tend to enlarge the extent of the necrotic region as the blood flow will evacuate the temperature inside the liver parenchyma, whereas if the tumor is closer to the vena cava, the advection will tend to reduce the size of the necrotic core as the temperature will be dissipated in the general blood flow vasculature. Therefore, we can claim that the effect of the advection is roughly proportional to the signed geodesic distance between the vena cava and the portal vein. Thus the impact of parenchyma perfusion is highly dependent on the tumor location.

Effect of the Probe Position

The position of the probe is a major unknown in our experiments. We assumed that the clinician put the probe at the center of the tumor (

Effect of the Segmentation

The segmentation also has an influence on the computed lesion (Table . 4.4 and illustrated on Figure 4.10). To notice the effect of user variability on the segmentation, the expert segmented twice the data of patient 03. This second segmentation was done independently and weeks apart from the first one. Figure 4.11 shows the superposition of the two segmentations. Based on these two segmented datasets, we computed twice the RFA ablation and compared the two computed lesions (Figure 4.12). It can be noted that the main difference between the two computations lies in the CFD results which is really sensitive to the segmentation, especially the segmentation of the vessels inlets.

Portal vein Tumor

Segmentation 1

Segmentation 2

Ground Truth

Computed lesions Figure 4.12: The two computed lesions for patient 03, the streamlines represent the parenchyma flow and are color-coded with respect to the velocity magnitude.

The total blood flow is the liver is the same (around 25 mL/s) but the parenchyma blood does not flow similarly in both cases due to the differences in the extremities segmentations.

Effect of Biophysical Parameters

The accuracy of the computation depends on the material properties since they are patient-specific, temperature-and space-dependent and they cannot be easily measured in vivo. Nominal tissue parameters were employed, with values often based on ex vivo experiments on animal tissue sometimes with a large varying range between published studies [START_REF] Tiesler | Identification of Material Parameters from Temperature Measurements in Radio Frequency Ablation[END_REF]]. A proper estimation of those parameters is needed but has often been overlooked due to its difficulty. A first step towards the personalization of those parameters is to perform a sensitivity study of the computation to the main parameters. From the results based on nominal parameters for patient 02, we can see that the heat conductivity is too low. For a heat conductivity 8 times greater than the nominal value ( 4.096 W (mK) -1 ), we manage to have a better outcome (Figure 4.14). Qualitatively, Figure 4.13 shows the variation in term of point-to-mesh error and similarity scores with respect to the diffusivity value. Quantitatively, Table 4.5 shows improvements in the point-to-mesh error (from 8.25 ± 5.92 mm to 3.85 ± 3.21 mm) and in the Dice index (from 61.7 % to 82.3 %).

Discussion

We have evaluated the first multi-physics model proposed for efficient patient-specific planning of RFA based on medical CT images. Despite possible biases in establishing correspondences between the post-to the pre-operative images due to registration or segmentation errors, and the use of nominal biological parameters, which are not patient-specific, our model provided promising results, opening new opportunities for RFA planning and guidance.

Effect of Large Vesssels

The computed lesions show that large vessels act as heat sinks as the lesion follows the large vessel's walls. A generic vena cava pressure, inflow and flow through the inlets of the hepatic veins were used for all computations. Since these values were not available, we assumed the pressure and the flow to be a constant in the physiological range of healthy patients. Yet, the pressure and the flow should be personalized as it can vary from patient to patient but also spatially and temporally inside the liver.

And the blood flow should instead be modeled as a pulsatile flux. These parameters could be adjusted based on catheter or imaging information such as Phase-Contrast MRI, which can measure the flow in the visible vessels. This study was performed with a weak coupling of heat propagation and CFD. Change in blood viscosity was neglected as we compute the blood flow only once, at the beginning and then we focused on the heat propagation and cellular necrosis only. A fully coupled computation where the blood viscosity and coagulation are considered could improve the outcome of the RFA model, i.e the extent of ablation and may provide insights on the entire physical mechanisms involved in RFA and tissue properties for long-term therapy prognosis. This aspect of RFA modeling will be investigated in future work.

Effect of Advection

The results of the fourteen ablations computed with the complete model and without the advection term demonstrate that the perfusion in the parenchyma does have an effect on the shape and the extent of the thermal induced lesion, depending on the location of the tumor. Therefore the perfusion in the liver may be an important factor to quantify the extent of ablation and predict therapies. Recent experimental studies on pig livers show an increase in the thermal conductivity with decreased distance to large blood vessels and with the perfusion rate and significant directional differences in thermal conductivity [Podhajsky 2009]. In this study, nominal values were used for the advection coefficient, which control the extent of the heat transport. Yet, this parameter could be estimated and animal data could help to better understand this phenomenon [Cheng 2002] using MRI-thermometry acquired intra-operatively for example.

Effect of the Probe Position

We showed that the probe position is an important unknown and the RFA computation could be improved by adding this information. As the probe insertion during the RF ablation is done under ultrasound guidance, 3D ultrasound images could be used to know the exact position of the probe. Additionally, the simulation of the probe using Dirichlet boundary condition is simplistic, and a geometrical model describing the exact shape of the probe would probably be beneficial for this framework.

Effect of the Segmentation

We saw that the advection and large vessels may have an important effect on the size, shape, and extent of the lesion. More extensive sensitivity analysis of this effect with respect to the segmentation of the vessels should be performed. Modification of the reconstructed anatomy, in particular in the circulation system tree may have an important impact on the parenchyma blood flow and therefore on the advection influence, as shown for patient 03 on Figure 4.12.

In this study, as only late venous phase CT images were available, we could not segment the hepatic arteries, thus the effects of the arterial flow were not considered to predict the extent of ablation. Yet, it is worth noting that the veins account for more than 70 % of the blood inflow of the liver parenchyma [Schenk Jr 1962]. A complete study would benefit from including the hepatic artery inflow, accounting for about one fifth of the hepatic blood inflow, but which comes at a higher pressure and pulsation characteristics compared to the portal vein inflow. If the three phases (arterial, portal and venous) CT images were available, it would be straightforward to include hepatic arteries as well as more venous vessels without any modifications for improved accuracy, our framework being modular.

Towards Personalization

Finally, tissue properties in the liver are spatially-varying and may also vary with the disease-state of the liver (e.g cirrhosis), but they are difficult to estimate from clinical data since they cannot be quantified from CT imaging for instance. In this study, all patients had normal liver parenchyma without cirrhosis or steatosis and heat propagation was modeled with standard tissue parameters reported in the literature for all patients. However our method makes it possible to estimate global patientspecific tissue parameters, which may increase the accuracy of the computation as it is the case in adjusting the heat conductivity for one patient in this study. Moreover, the perfusion of the tumor was not taken into account which can impact its thermal conductivity. Yet, the border zone beyond the margin of the tumor is of greater interest and is considered as regular hepatic tissue. Further investigations and experiments on patient-specific tissue parameters will be necessary to include these new findings into our model.

Preliminary experiments suggest that systematic sensitivity analysis and personalized computations with patient-specific and not nominal parameters are necessary to get a better match with the actual outcome [Audigier 2014a].

Conclusion

Before deploying model-based therapy planning or guidance in clinical settings, clinical validation is a necessary step. The computational model of RFA needs to be validated on extensive data from a large population dataset. The proposed framework make it possible due to its level of integration. Our system may thus constitute a first step towards clinical application of a RFA computational model.

From a research point of view, the proposed framework paves the way to the quantitative and systematic evaluation of computational models of RFA. There is growing evidence that patient-specific anatomical models are necessary to fully comprehend RFA treatments, in particular in the light of recent in vivo studies [El-Serag 2003]. By providing a complete system for patient-specific modeling, our system may be used as input to more comprehensive inverse problem studies, and constitute a useful surrogate tool for RFA planning, potentially improving the outcomes for the patients.

Introduction

Motivated by the promising results obtained by fitting the diffusivity in the last chapter (section 4.4.5), we present in this chapter, a framework for the personalization of RFA which enables patient-specific parameter estimation in the case of multiple hepatic lesions. The forward model relies on the computational model of heat transfer and cellular death during RFA, presented in Chapter 3 and it is based on anatomies estimated from CT (Section. 3.2.3). The computation model of the blood flow is the Fully Coupled LBM CFD model so the Lattice Boltzmann Method (LBM) is used to compute not only the heat diffusion and cellular necrosis as in [Audigier 2013a] but also blood and parenchyma flow in the liver tissue. This framework is particularly efficient for the personalization as it provides a fast solver and naturally accounts for the flow transition between veins and parenchyma. The effect of the small vasculature inside the liver which are not visible on the CT data is taking into account.

The model is then personalized based on a first ablation. Our framework is adapted to situation where no temperature map is available. This information is used to plan subsequent ablation(s) of the same or additional lesions to treat for a more refined, targeted ablation that minimizes probabilities of recurrence. This can be validated in case of the ablation of multiple tumors inside the liver, assuming that the parameters are spatially invariant within the same patient.

In Section. 5.4, heat conductivity and porosity were selected as the most sensitive parameters for predicting the necrosis extent. After their estimation on patientspecific data, we demonstrate improved prediction accuracy. We discuss the contributions of this chapter and perspectives in Section. 5.5.

Available Clinical Data

In this chapter, the personalization method used constrains us to consider only patients who had several tumors ablated. Our database is a subset of the database presented in Section. 4.2, which comes from the Johns Hopkins Hospital in Baltimore, MD, USA. Our database is made of CT scans from three patients who underwent RFA of several hepatic tumors. For each patient, pre-and post-operative late venous phase CT are available. We generate volumetric binary images of the parenchyma, tumors, hepatic veins, vena cava, portal vein, without the hepatic artery since only single-phase CT images are available.

Mathematical Model of RFA Simulation

Forward Model

The forward model used in this chapter relies on the computational model of RFA, presented in Chapter 3. The simulation of heat transfer and cellular death inside the liver depends on the patient-specific anatomy (Section. 3.2.3). The nominal formulation from [ONeill 2011] of the cell death model (with α = 1) is used to compute necrotic extent. The computation model of the blood flow is the Fully Coupled LBM CFD model so the Lattice Boltzmann Method (LBM) is used to compute not only the heat diffusion and cellular necrosis but also the blood flow inside the main vessels and the parenchyma considered as a porous medium (Section 3.2.2). This latter method is based on a computational fluid Dynamics (CFD) solver which incorporates a porous part to deal with the liver parenchyma. The different steps of our method are illustrated in Figure 5.1. The main parameters of the model of heat transfer and cellular necrosis in liver tissue are optimized to match the observed extent of the necrosis.

Parameter Estimation

As parameters in the heat transfer and cellular necrosis equation are customarily taken from the literature, we aim to personalize them given the observed extent of the necrotic region measured post-operatively as temperature maps are not readily available. Most of the parameters are defined as constant whereas the heat capacity c t and conductivity dt = dt * (1 + 0.00161) * (T -310)W (mK) -1 are temperature dependent and therefore spatially distributed [Audigier 2013a]. Using DAKOTA 1 , we first perform a sensitivity analysis to know which parameters mostly influence the volume and the point-to-mesh error [START_REF] Zheng | [END_REF]] of the computed necrosis area. Then, we optimize the most sensitive ones: the heat conductivity and the porosity, as to minimize the average point-to-mesh error between computed and observed necrotic region. To this end, we use a gradient-free optimization method, the Constrained Optimization BY Linear Approximations (COBYLA), which required only a few numbers of forward simulations.

Global Sensitivity Analysis: Variance-Based Decomposition

In order to know which parameters of the model mostly influence the output of the RFA simulation, we first perform a global sensitivity analysis using a variance-based decomposition.

Variance-based Decomposition (VBD) is a global sensitivity method that summarizes how the uncertainty in the model output can be apportioned to uncertainty in individual input variables [START_REF] Adams | [END_REF]]. VBD uses two primary measures called the Sobol indices: the main effect sensitivity indice S i and the total effect indice T i .

Let us call Y, the given output of the computation (in our case, it is either the volume of the necrosis or the point-to-mesh error between the simulated necrosis and the ground truth). Y is a function of m input variables: Y = f (x). We also consider: x -i = (x 1 , ..., x i-1 ; x i+1 , ..., x m ). The Sobol indices can be computed as follow:

Si = V ar xi [E(Y |x i )] V ar(Y ) (5.1)
and:

T i = E(V ar(Y |x -i )) V ar(Y ) = V ar(Y ) -V ar(E[Y |x -i ]) V ar(Y ) (5.
2)

The main effect sensitivity index, S i corresponds to the fraction of the uncertainty in the output, Y , that can be attributed to input x i alone. The main effect sensitivity index compares the variance of the conditional expectation V ar x i [E(Y |x i )] against the total variance V ar(Y ). Whilst, the total effect index, T i corresponds to the fraction of the uncertainty in the output, Y , that can be attributed to input x i and its interactions with other variables. It should be noted that: m i=1 T i ≥ 1, due to the fact that the interaction effect between two different input variables is counted twice [START_REF] Adams | [END_REF]].

Experiments and Results

All experiments were executed on a Windows 7 desktop machine (Intel Xeon, 2.80 GHz, 45GB RAM, 24 CPUs) with a Nvidia Quadro 6000 1.7 GB (448 CUDA cores). A synthetic case, presented on Figure 5.2, is used to validate the CFD solver. A cylindrical mesh with spherical inlet and outlet is used. Blood flow is set at the inlet and pressure at the outlet. Two cases are considered. The first case presented on (Figure 5.2, Left), simulate a single vessel, with porosity set to 1 inside the vessel, 0.04 in the endothelium and 0.1 outside. In the second case, presented on (Figure 5.2, Right), the porosity map is used to emulate two veins and the liver parenchyma. The mass conservation is verified by comparing the flow through the surface of the "vessel" at z = 20mm: Φ in and at z = 80mm: Φ out . The computation is run in both cases with a spatial resolution of 0.5 mm and until a steady state is reached (tolerance: 10 -4 ). In the first case, the difference between these two fluxes is 3%, whereas in the second case, 2.8% difference between these two flux is achieved.

Evaluation of the CFD Solver on a Synthetic Case

Sensitivity Analysis

We want to know the sensitivity of four uncertain parameters of our model: dt, H, kf , ε on the volume of the computed necrosis but also on the point-to-mesh error between the computed necrosis area and the one computed with the nominal parameters from the literature. To this end, a synthetic case has been setup to speed-up the process (Figure 5.3). The range of parameters values used [Altrogge 2012] are reported in Table 5.1. These parameters of interest were modeled with a uniformly distributed uncertainty, and the sensitivity analysis was performed using variance based decomposition to compute the global sensitivity indices (so-called Sobol indices). dt has the largest total effect (0.58) as compared with kf , H and ε (0.16, 0.15, 0.43) on the volume of the lesion, whereas ε and dt have the same larger total effect (0.37) as compared with kf and H (0.16, 0.35) on the point-to-mesh error with respect to Notation Parameter Name Min -Max H convective transfer coefficient 24.4×10 4 48.8×10 4 W (m 3 K) -1 dt heat conductivity 0.25 1.24 W(m K) -1 kf damage rate coefficient 3.2×10 -3 3.4×10 -3 s -1 ε porosity 0.1 0.9

Verification of the Optimization Framework

In order to confirm the accuracy of the optimization framework, we considered a synthetic case on a regular cuboid domain (Figure 5.3) where all the phenomena occuring during RFA were present (diffusion, reaction and advection) and with nominal parameters of tissue properties. First we simulated a necrotic area with generic parameters by emulating the clinical RFA protocol: during 7 minutes, we heated at 105 • C, the simulation continued for 3 more minutes without heating so that each cell reached a steady state. Then, the main parameters: dt and ε are estimated by minimizing the mean of the point-to-mesh error between the computed necrotic region and the one created at the first step. We managed to obtain the estimated parameters with 6.1% of error on dt and 2.1% on ε in 32 minutes after 36 iterations with a mean of the point-to-mesh error of 10 -3 mm. The convergence of the COBYLA algorithm is illustrated on Figure 5.4.

Evaluation on Patient Data

We evaluated our model on 3 patients, with 7 ablations (several tumors ablated for each patient) for whom pre-and post-operative CT images were available. Clinical RFA protocol was simulated: the probe was deployed within the tumor and cells in a diameter defined pre-operatively around the center of the tumor probe tip were heated at 105 • C during 7 minutes or 2 times 7 minutes. The diameter and heat duration were iteratively increased according to the size of the tumor. The simulation continued for 3 more minutes without any heat source so that each cell reached a steady state. The parameters were considered spatially invariant within the same patient. To the best of our knowledge, all the proposed computational models of RFA assume numerous simplifications, such as the homogeneity and isotropy of the tissue. All the characteristics (electrical, thermal, etc) are normally considered to be isotropic. Nevertheless, the cell death model locally changes the properties of the tissue, and different parameters are related to different location inside the liver.

For instance the heat capacity of the tissue (c t ) is different from the heat capacity of the blood (c b ). Moreover, the way the heat is propagating depends on the location inside the liver. There is a reaction term (H, convective transfer coefficient) in the large vessels, and an advection term in the parenchyma, which depends on the porosity ε, but we consider that they have a constant value. The parameters are estimated on one tumor by reducing the error with the ground-truth. Then, we computed the cell death area of the other tumor(s) of this patient with the personalized parameters. The computed cell death area compared qualitatively well with the observed post-operative necrosis zone for tumor located at different place inside the liver, close to large vessels, or on the border. The predicted lesion was manually segmented by an expert and registered to pre-operative image. Quantitatively, average point-to-mesh errors (Table 5.3) were within tolerance in clinical routine for the four tumors estimated with the personalization of the main biological parameters, as the probes can be deployed in steps of 1 cm. Figure 5.6 shows the error in term of point-to-mesh errors and Dice scores for the 7 ablations simulated with personalized parameters. The estimated heat conductivities were lower than the nominal value (0.512 W(m K) -1 ), whereas the porosity was very close (0.1), as illustrated on Figure 5.7. Patient 2 presents a quite large Dice difference between the two cases. It might be due to segmentation or registration issues, and potentially to model limitations (assumptions made, etc.), which will be further investigated in a pre-clinical setup.

The main parameters are first considered with nominal values taken from the literature and then optimized to match the observed extent of the necrosis. Current errors can be explained by segmentation and registration processes but also by the limited number (2) of personalized parameters.

Other experiments on patient 1 showed also a significant improvement of the correlation between predicted and actual ablation extent compared to the prediction using only nominal parameters (average point-to-mesh errors of 4.44 mm vs 4.98 mm, average Dice score of 72.0% vs 68.5%). Those results are illustrated on Figure 5.8, Figure 5.9 and Figure 5.10. The comparative results with simulations with nominal parameters for patient 1 are reported on Table 5.4 and Table 5.5.

Ground truth

Vena

Discussion

The personalization of the sensitive tissue parameters allows to have a better estimation of the necrosis and to predict the outcome of RFA in case of multiple tumors inside the liver. As we are using the Fully Coupled LBM CFD solver, our framework totally rely on LBM, no advanced meshing techniques are required. All the computations are directly done from patient images: heat propagation and cell death modeling as well as the heat sink effect of blood vessels and porous circulation in the liver. The current method needs several tumors for validation and is worth using only when no temperature maps are available, but it could be easily translated into clinical settings. Adaptation for RFA under image-guidance is considered: RFA procedure is usually done in several steps (increase in probe diameters for example). An intra-operative image is acquired at the end of the first step and used to personalize key parameters, providing a powerful guidance tool. No post-operative images are required. A necessary step before deploying this method in clinical settings is a pre-clinical validation with extensive data on larger populations to evaluate the

Introduction

In this chapter, we introduce a pre-clinical study for the validation of RFA model, based on pre-, intra-and post-operative data.

The sophisticated experimental protocol is described first (Section 6.3), followed by the RFA model we aim at validating. The computational model relies on subjectspecific anatomies estimated from Computed Tomography (CT) (Section 6.4). Heat transfer is computed according to the bio-heat equation (Eq. 2.4) and the Fully Coupled LBM CFD Model is used: a Computational Fluid Dynamics (CFD) solver incorporates a porous part to deal with large vessels and the liver parenchyma hepatic blood flow simultaneously as described in Section 3.2.2. Moreover, this is coupled with a novel cellular necrosis model.

Finally, we present the data analysis required before using the information from the pre-clinical images into the model (Section 6.5). Section 6.6 concludes this chapter.

Available Pre-clinical Data

Our database is made of CT scans, MR images, invasive and device-based measurements from five pigs who underwent RFA of one or several hepatic surrogate tumors. Figure 6.1 illustrates the case of Pig 5. For each pig, pre-, intra-, and post-operative images are available. We generate volumetric binary images of the parenchyma, tumors, hepatic veins, vena cava, portal vein and the hepatic artery since three phases CT images are available. The pre-clinical data were processed off-line. The results were not available to the interventionist during the RFA procedure. This set of swine data has been undertaken at the IHU of Strasbourg in France.

Surrogate Tumor

Post-op lesion

Gall Bladder 6.3 Pre-clinical Experiment

A Pre-clinical Study for RFA Model Validation

A comprehensive animal experiment has been realized for model validation and parameter identification including several modalities at pre-, intra-and postoperative stages. Figure 6.2 illustrates the available pre-clinical data (in blue), the different pre-processing steps needed (in red, green and white) before performing the computation (in grey). We have developed this pre-clinical study specifically for the validation of the model developed in chapter 3. We have seen that this model is based on patient images. For this reason a pre-operative CT scan of the liver with the implanted surrogate tumors is acquired in order to get an accurate representation of the anatomy. We showed in chapter 4 that the large vessels play an important role in the shape of the lesion, so we acquired pre-operative CT images with the three available phases including portal, venous and arterial phases. One limitation of the proposed model is the probe representation and its location, as described in section 4.5.3. To tackle this issue, a CT scan of the probe alone was first acquired to get an accurate model of the RFA probe, and then an anatomical MR image was acquired after the insertion of the probe to know exactly its position. The boundary conditions of the CFD solver were fixed with nominal flow and pressure values coming from the literature, as described in section 4.5.1. To get a more personalized computation, Phase-Contrast MR images were acquired to get a subject-specific blood flow and pressure measurements from catheter were performed to get invasively the subject-specific pressure. The clinical protocol employed for each specific patient was unknown, so in this study, we decided to register for each ablation procedure the delivered power and the temperature distribution at the thermistors provided by five tips of the probe. Finally, as the evaluation of the model relies on the comparison with the actual necrotic lesion segmented on the post-operative image, we decided to acquire post-operative images of the lesion two days after the intervention.

Experimental Set-up

The present experimental study (No. 38.2014.01.063) received approval from the local Ethical Committee on Animal Experimentation. All animals used in the experimental laboratory were managed according to French laws for animal use and care and according to the directives of the European Community Council (2010/63/EU).

Pigs are considered as a relevant animal model as their hepatic system is similar to the human one. Several surrogate tumors (diameters < 3cm) are inserted on five swines at various locations of the liver (close to vessels or the Glisson capsule) under ultrasound (US) guidance (SIEMENS ACUSON S3000), followed by the acquisition of pre-operative CT images including portal, venous and arterial phases (CT scanner: SIEMENS Somatom Definition A5). The surrogate tumors are made of a specific gel (a mix of biocompatible gelatin, alginate and nanoparticles) which exhibits an hyper-intense signal in CT and MRI as illustrated in Figure 6. The surrogate tumors are mainly made of gelatin, whose thermal properties can be found in the literature. In [START_REF] Culjat | [END_REF]], the gelatin density is reported to be 1050 kg.m -3 , which is close to the value used in this work (1060 kg.m -3 for blood and tissue density). In [Dunmire 2013], the gelatin heat capacity is reported to be 3.6 kJ(kg K) -1 , which is exactly the value used for the heat capacity of healthy tissue in this work. However, the reported value of conductivity is approximately 0.303 + 0.0012 (T -310) W(mK) -1 if T is in Kelvin , which differs from the value used: 0.512 (1 + 0.00161 (T -310)) W(mK) -1 . In order to know how the surrogate nature of the "tumors" impacts the heat transfer and subsequent "necrosis", it could be possible to run the experiments by setting this specific value at the points belonging to the surrogate tumor. However as we used the conductivity value of the tissue for the tumor cells in the clinical study (Chapter 4), similarly we kept this assumption in the pre-clinical study.

Invasive catheter is introduced through the jugular vein to get the free and the wedge pressure in a subhepatic vein and in the vena cava. An MR-compatible RFA probe, the radiofrequency interstitial tumor ablation (RITA) probe (StarBurst RFA, AngioDynamics; www.angiodynamics.com), is deployed at 2 cm of diameter (the diameter of the area defined by the tips of the probe is 2cm) under US guidance (Figure 6.4) next to the targeted surrogate tumor. An MR image is then acquired to get the position of the probe in the liver (SIEMENS Magnetom Aera 1.5T) and flow data. The temperature and delivered power are monitored and recorded intraoperatively during and after the ablation. Finally, a post-operative CT, T2 or T1 + gadolinium MRI is acquired two days after the ablation to assess the extent of the necrotic areas. Overall, pre-, intra-and post-operative images are available, along with interventional device measurements (Figure 6.2). To the best of our knowledge, no such validation set-up has been reported previously in the literature.

RFA probe liver

RFA probe tumor

Computational Experiment

Overview

Model validation can be performed, as it is done in chapter 4 through a comparison between the computed necrotic extent and the actual lesion observed on the postoperative image. We will proceed as in chapter 4, where all the computations are performed on the pre-operative image, and the accuracy of the computation is evaluating by a comparison with the post-operative image. We rely on the RFA computation described in chapter 3 to compute the heat transfer from the bio-heat equations (Eq. 2.4). To compute the hepatic blood flow, the Fully Coupled LBM CFD Model described in Section 3.2.2 is used as we want to obtain simultaneously the blood flow in the large vessels, in the arteries and in the parenchyma. Chapter 5 also showed the importance of patient-specific parameters. In this chapter, we rely on intra-operative data (the delivered power and the probe temperature) to show that subject-specific parameters can be estimated.

RFA Modeling

Blood Flow Computation

The Fully Coupled LBM CFD Model described in Section 3.2.2 is used to compute simultaneously the flow in the main vessels, in the arteries and in the parenchyma. This model depends on the porosity coefficient ε (fraction of blood volume over the total volume) defined in this chapter through a porosity map (Figure 6.13) as described later in the next section (section 6.5.3.4).

Boundary Conditions

This method makes the boundary conditions simple to treat: no boundary conditions are fixed on the extremities of the vessels inside the parenchyma thanks to the use of the porosity map, contrary to [Audigier 2015a]. At the border of the liver, no flux boundary conditions are used (Neumann) whereas Dirichlet boundary conditions are applied at the inlets of portal vein, vena cava and hepatic artery and at the outlet of the vena cava and hepatic artery. Pre-operative Phase-Contrast MR images are used to set the portal vein, vena cava and arterial inflows, ϕ p , ϕ vc in and ϕ a , whereas invasive pressure measurements through catheter are used to set the vena cava and arterial outlet pressure p 0 and p a . Figure 6.5 illustrates boundary conditions on a subject-specific geometry.

Heat Transfer Model

In this chapter, the heat transfer model described in chapter 3 is used. It depends on the patient-specific anatomy (the segmentation process is presented in Section 6.5.2) and on the blood flow inside the main vessels and the parenchyma considered as a porous medium (Section 6.5.3.4). The Pennes Model is solved in the large visible vessels and the WK model in the liver parenchyma. Nominal values of parameters reported on Table 6.1, are used.

Adapted Formulation of Cellular Necrosis Model

The three-state model described in section 3.2.2.2 is used to compute tissue necrosis. Unlike in [ONeill 2011], the constant α is different than 1 ( α = δ/ β(1 + 10N )) to

.3×10 -3 γ (s -1 ) recovery rate coefficient 7.7×10 -3 δ (s -1 ) vulnerable rate coefficient 3.3×10 -3 c V t (J(kg K) -1 ) heat capacity of vulnerable tissue 3.6 ×10 3 c N t (J(kg K) -1 ) heat capacity of necrotic tissue 0.67×10 3 c U t (J(kg K) -1 ) heat capacity of undamaged tissue 3.6 ×10 3 dt (W(m K) -1 ) heat conductivity 0.512 R (W (m 3 K) -1 ) reaction coefficient 27.1×10 4 α v advection coefficient 0.11
decouple the damage rate coefficient from the vulnerable rate coefficient, so that three distinct transition rates are considered to allow cells to reach the vulnerable state. Figure 6.6 represents the solution of Eq. 3.11 at one vertex over time if a constant temperature of 105 • C is applied.

Numerical Computation using LBM

For a time step of ∆t = 0.5 s and a spatial resolution of ∆x = 0.5 mm, better than real-time computation can be achieved on a desktop machine (Windows 7, Intel Xeon, 3.30GHz, 16GB RAM, 12 CPUs, Nvidia Quadro K5000 4.0 GB). For example, it took 577 seconds to compute 1000 s of ablation of tumor 4-1.

Parameter Estimation from Probe Measurements

During the intervention, the delivered power and the temperature distribution are measured by the ablation probe itself. We explore how this can be used to estimate the parameters of the model. During the heating phase, the simulated heat power P s can be compared to the measured one P m . During the cooling phase, the simulated temperatures T s can be compared with the measurements T m , read from five tips of the probe (four tips do not have any thermistors).

In this study, the parameters in the heat transfer and cellular necrosis equation come from the literature. After a sensitivity analysis reported in the next section 7.4, we choose to estimate the heat capacity c U t and the constant part of the conductivity dt as they mainly influence the delivered power, the temperature distribution and the size of ablated regions.

As temperature maps are not readily available, these two parameters are personalized from probe measurements (temperature and power) by minimizing the following cost function:

arg min c U t , dt 1 t a t=ta-dtm t=0 (P m (t) -P s (t)) 2 σ 2 Pm + 1 t c -t a t=ta+tc t=ta (m Tm (t) -m Ts (t)) 2 σ 2 Tm (6.1)
where dt m is the time step of the measurements, σ P m and σ T m are the standard deviations associated with the heat power and the temperature, both of them evaluated from the variability in the available observations (equal to 13.3 W and 5.1 • C in our experiments). σ P m is the standard deviation of the power of all the twelve ablations at one time point during the heating phase (t= 50 s), which acts as a normalization factor in Eq. 6.1. We chose t = 50 s to define σ P m since all the power curves are coherent at that time. Moreover, as illustrated in Figure 6.7, this time is a good compromise since it is not too early, avoiding the boundary conditions effect, and not too late, before the cell death could affect the delivered power. Eight ablations are performed through several short cooling and heating periods, in this case, the delivered power is reduced during the different cooling stages. On the other hand, the other four ablations included only one long final cooling stage after a continuous heating period (red, blue, cyan and green curves in Figure. 6.7). Figure . 6.8 shows the four cases whose the delivered power over ablation is a long heating phase followed by a cooling stage, and one typical case whose ablation is performed through several short cooling and heating periods (pig 5 -tumor 1). σ T m is the standard deviation of the temperature computed from one ablation during the cooling phase (we did not use all the twelve ablations as their cooling stage did not start at the same time), where all the temperature curves at the 5 tips of the probe are coherent. It is defined as follow:

m T (t) = 1 5 5 i=1 T i (t) σ T 2 m = 1 t c t=ta+tc t=ta ( 1 5 5 i=1 (T i (t) -m T (t)) 2 )
To cope with the uncertainty in the rotation of the probe along its axis, the mean tip temperature m T (t) is used for the personalization instead of directly mapping the tip temperatures.

Data Pre-processing

The model is based on a representation of the anatomy obtained from pre-operative images. It takes as inputs the geometry and the location of the probe and as boundary conditions blood flows and pressures. Finally the validation of the model outcome is based on a comparison with the actual lesion observed post-operatively. In order to get this information, data pre-processing is required. Segmentations of the RFA probe and of subject-specific anatomy coming from the pre-and postoperative images are needed. The position of the RFA probe is known from an intraoperative anatomical MR image, but its registration to the pre-operative image is not straightforward, given the artifacts induced by the RFA probe during the acquisition. Moreover, the segmentations of the visible vessels come from anisotropic images. To get stable solution of the CFD solver, smooth vessel trees have to be generated. Since the acquisition process of the Phase-Contrast MR images was complex, preprocessing of the images was necessary to get coherent settings of the boundary conditions. Finally the registration of the post-operative image to the pre-operative image was extremely challenging. It required the use of a robust and advanced registration algorithm. We present here one registration method which gave the most significant results out of several methods which have been tried.

3D Modeling of the RFA Probe

The position and the geometry of the probe are inputs of the RFA model, thus we need to extract these information from images. The geometry of the probe deployed at 2 cm is acquired from a CT image (resolution: 0.2x0.2x0.9 mm) of the probe alone. A 3D mesh is then reconstructed by thresholding the image intensity (Figure 6.9, left panel), and manually registered to the pre-operative CT using the main axis of the probe and intra-operative MR data. The manual rigid registration is done using Paraview1 and is visually checked by an expert up to a rotation along the probe axis while the main axis is clearly visible as well as three or four tips (depending on the ablation considered) but the MR resolution does not allow to distinguish between the nine tips of the probe individually (Figure 6.9, right panel).

Geometrical model of the RFA probe

Tips of the probe

Segmentation

As each structure of interest is needed as input of the RFA computational model, the segmentation of patient images is a determinant task. From the pre-operative CT data, the following anatomical and pathological regions are segmented semiautomatically and meshed by the tools of Visible Patient, Strasbourg, France (Figure 6.10): parenchyma, hepatic veins, vena cava, portal vein, hepatic arteries and all tumors. These regions are then used to define the computational domain. The meshes are rasterized and a multi-label mask image is created to identify the different structures. To define the computational domain, a level set representation of the liver without tumor and vessels is computed. From post-operative CT or MR data, necrotic areas are segmented and meshed as well.

CFD Computational Domain

Vessel Trees

The pre-operative images from which the vessels are segmented have anisotropic resolution. To avoid unstable solution of the CFD solver, smooth vessel trees have to be generated from the semi-automatic segmentations. To this end, centerlines are extracted from each vessel segmentation using VMTK2 (Figure 6.11, left panel). Along each branch of the centerline, the mean radius is computed and the smooth vessel tree is generated by overlaying cylinders having this mean radius in a piecewise fashion (Figure 6.11, right panel). All those cylinders are rasterized on a single image.

Portal vein

Vena Cava

Centerlines

Smooth vessels: the color encodes the mean radius of each tube 2D+t Phase-Contrast MRI have been acquired pre-operatively before the probe implantation. For each pig, several 2D+time sequences are acquired at different time points, for the reproducibility of data, but also at different locations: at the inlet of the hepatic artery, at the inlet of the portal vein and at the inlet of the vena cava. Given one 2D+t Phase-Contrast MRI sequence acquired at one time point, the user places a single seed in the vessel of interest in the first image and then an automatic method is used to segment this vessel on each 2D slice (Figure 6.12, right panel) [START_REF] Gulsun | [END_REF], Guetter 2011]. Using the vessel area defined by the segmentation, the mean blood flow can be computed at each time. If N 2D+t sequences are acquired at the same location but at N different time points, this is done N times and the RMS of the mean of those N curves is used. The acquisition process of flow is complex to handle and requires a significant learning curve. Due to the nonreproducible and non-coherent measures on the first four pigs, the measures of Pig 5 are used for the five pigs in all the computations (Figure 6.12). We assume that flow would not vary so much from pig to pig. It is a reasonable assumption, since it was demonstrated that liver weight and hepatic blood flow in all mammalian species could be readily related to body weight by a simple equation [START_REF] Boxenbaum | Interspecies variation in liver weight, hepatic blood flow, and antipyrine intrinsic clearance: extrapolation of data to benzodiazepines and phenytoin[END_REF]]. The 5 pigs used in the experiments were all healthy, of similar age (3-4 months) and weight (33, 29.6, 26.7, 25, and 28.4 kg), that is why we assumed that the blood flow would not vary so much from pig to pig.

For Pig 5 (vena cava), the acquisition has been done at 3 different time points (N = 3). The RMS of the mean of those three curves is used (Figure 6.12, left panel) since numerical results have shown that pulsatile velocity profile in large blood vessels has little difference in effect on the thermal lesion region of tissue compared with uniform or parabolic velocity profile [Horng 2007]. Figure 6.12: (Left): Velocity in the vena cava of Pig 5 from Phase-Contrast MRI with respect to trigger time (all measures are registered on the cardiac pulse). In cyan, the mean of different measurements at different times (red, blue, green curves). The RMS of this curve (black curve) is then used as boundary conditions in the computation as ϕ vcin . (Right): Phase-Contrast MRI of Pig 5: the vena cava is segmented in green, the portal vein in blue and the aorta in yellow.

Pressure from invasive measurements

The computational model of RFA requires also vena cava and hepatic artery blood pressures as inputs. In order to impose subject-specific values as boundary conditions, these blood pressures have been measured to avoid the use of nominal values from the literature.

The pressures at the outlet of the vena cava and the hepatic artery have been measured invasively by catheter introduced through the jugular vein. The same values (from Pig 5) are used as boundary conditions in the five pigs, as it is the case for the blood flow measurements.

Porosity Map

The porosity has to be defined everywhere in the computational domain as it is an input of the CFD solver. The vessel walls are defined as follow. The smooth vessel trees are rasterized on a single image. On this image, we perform a 26-connexity dilation on the voxels of the vessels to model the endothelium and avoid the flow to go through the vessel wall. The extremities are detected using the centerlines previously used and the porosity at the extremities is set to the parenchyma porosity value. Figure 6.13 shows an example of porosity map. The porosity value is 1 in the CT-visible vessels, 0.1 in the porous parenchyma [Payne 2011], modeled then as a porous medium. Experiments are performed to obtain a sufficiently small porosity (0.04) in the vessel walls (impermeable medium) to avoid the occurrence of shear stress on the vessel walls (leakage). Figure 6.13: (Left): The porosity map of Pig 1 created from the vessels segmentation. The porosity is 1 in the vessels (in white), 0.04 in the endothelium (in black) and 0.1 in the parenchyma. (Right): Zoom on the porosity map. The extremities of the vessels are detected using the centerlines so that the blood flow can go through the vessel extremities.

Post-To-Pre Registration

Due to ethical reasons, CT with contrast agent could not be acquired two days after intervention, making accurate post-to pre-operative registration challenging. In order to compare the results of the computation with the ground truth given by the post-operative images, registration of the post-operative image to the pre-operative image is performed. For each pig, the pre-operative image (CT from different phases) does not necessarily belong to the same modality as the post-operative image (CT or MRI (T2 or T1 + gad)) due to individual experimental changes. Moreover, the pigs are neither in a similar position nor a similar condition, as they had an empty stomach the day of the intervention, which was not true two days after. Breathing has also an important impact on the shape of the liver. For those reasons, we choose to register the post-operative meshes to the pre-operative image instead of pure image-to-image registration. From the semi-automatic segmentation of the post-operative image, a tetrahedral mesh is generated using CGAL3 . From the segmentation of the pre-operative image, a binary mask of the liver including vessels and arteries is generated. Using the SOFA framework4 , the post-operative mesh is deformed elastically to match the intensity profile in the pre-operative image using the Finite Element Method (FEM) and a co-rotational model [Faure 2012]. First-order implicit Euler time integration is employed and the system of equation is solved with the conjugate gradient algorithm. Figure 6.14 shows the preoperative binary image with the surface of the post-operative mesh before the registration (Figure 6.14, Up) and after the registration (Figure 6.14, Down).

Before Registration After Registration

The computed deformation field is then applied to the necrotic surface mesh using a barycentric mapping between the coordinates of the surface mesh and the tetrahedral mesh. The registration is fast to perform. As an example, it took 110 seconds for Pig 5 to register the post-operative mesh to the pre-operative image.

Conclusion

In this chapter, we presented a sophisticated pre-clinical set-up to validate a complete and complex multi-physics model of radiofrequency ablation (RFA). This comprehensive validation is based on pre-, intra-, post-operative images and device-based measurements. To our knowledge, it is the most complete validation set-up proposed in the literature. The RFA computational model considered takes into account the main occurring biophysical phenomena (heat transfer, cellular necrosis, hepatic blood flow). The advection effect of the porous circulation in the parenchyma [Audigier 2013a] is taken into account by including a CFD model, robust to image noise and anisotropy to compute the venous and arterial blood flow in the liver parenchyma. By adding such a component, the proposed model is further improved.

The pre-clinical validation approach is difficult to establish. The experimental set-up requires different acquisitions of different modalities (CT, MRI, US), at different times (pre-, intra-, post-operative images). The different images at Day 1 (day of intervention) should be acquired quickly, as the anesthesia of the pig cannot be too long since the pig will be kept alive for 2 additional days. Due to the high complexity of the experimental set-up, pre-processing is necessary. The segmentation was performed semi-automatically, but the vessel meshes had to be smoothed in order to avoid unstable solution of the CFD solver. The probe was segmented from a CT image with a good resolution, the probe geometry considered in the model is convincing. However the artifacts induced by the probe and the bad resolution of the interventional MR images, did not allow to accurately register the probe position. The setting of subject-specific boundary conditions for the blood flow and the pressure was not straightforward, as the acquisition process of Phase-Contrast MR images was complex to handle and required a significant learning curve. Due to the non-reproducible and non-coherent measures on the first four pigs, the measures of blood flows and pressures of Pig 5 only were used for the five pigs in all the computations. Finally, after several trials of registration algorithms, we were able to reliably estimate the deformation between pre-and post-operative imaging for a precise validation of the necrosis extent prediction. The validation of the RFA model on the pre-clinical data presented in this chapter is performed in the following chapter.

Introduction

In this chapter, the pre-clinical study presented in chapter 6 is used for the evaluation of RFA model. Since our approach relies on pre-clinical data of healthy pigs in a controlled environment, the evaluation leads to an increased confidence in the computed information.

The computational model evaluated here was presented in the last chapter (Section 6.4). An evaluation of each pre-processing steps is performed first (Section 7.2). In Section 7.6, the model is evaluated on twelve ablations from five healthy pigs, to which surrogate tumors have been implanted. We demonstrate accurate predictions on pre-clinical data both for temperature, power evolution and necrotic region extension: the mean error between measured and simulated temperature (respectively power) is 5.1 • C (respectively 25.9 W), which is good compared to the targeted ablation temperature of 105 • C and the maximal power of 150 W. The mean pointto-mesh error between predicted and actual ablation extent is 5.3 ± 3.6 mm, which is smaller than the diameter of the targeted surrogate lesions (around 2 to 3 cm). As we have a LBM implementation on graphics processing units (GPU), a single simulation is faster than real time, which allows us to adjust model parameters for a better evaluation in Section 7.6.2. Key biophysical parameters (heat conductivity and heat capacity) are estimated from this pipeline to show that parameter personalization could be achieved. We use inverse modeling to minimize the error between the computed power and cooling temperature and the observed values (the delivered power directly measured by the RFA device itself and the observed temperature drop during cooling measured at the 5 tips of the probe). Section 7.9 concludes this chapter.

Evaluation of the Pre-processing

Evaluation of the Registration

First, in each case, the accuracy of the registration was visually checked by an expert as it is illustrated on Figure 7.1. Figure 7.1: Results of the Post-to-Pre registration on the liver of Pig 1. (Up): Superposition of the post-operative liver in red with the pre-operative liver before the registration; (Down): Superposition of the post-operative liver in red with the pre-operative liver after the registration. Table 7.1: Evaluation of the registration: Point-to-mesh errors from the preoperative meshes to the post-operative meshes in millimeter before and after the registration for Pig 2 are computed. Before the non-rigid registration, all the meshes are centered.

Mesh

Before Whenever available, the non-rigid transformation between pre-and postoperative images is applied to the post-operative vessel and/or tumor meshes and we can compare them to their pre-operative equivalent meshes. Point-to-mesh errors from pre-operative meshes to registered post-operative meshes are computed and reported in Table 7.1 for Pig 2. Before the non-rigid registration, all the meshes are centered. The error is decreased after the non-rigid registration except for tumor 1 where the two meshes are already in good agreement before the registration. As illustrated in Figure 7.2, the point-to-mesh errors are high for the arteries as the segmented vessels differ in both segmentations of Pig 2.

Mesh from

Pre-op image Before registration After registration Figure 7.2: Comparison between the pre-operative arterial mesh in blue and the post-operative meshes (in red: before the registration and in green: after the registration) of Pig 2. It is clear that the arterial vessels segmented on the pre-operative mesh differ from the one segmented on the post-operative mesh. However, the main part of the vessel (vertical part) is correctly registered. Figure 7.3 shows the results of the registration algorithm on the liver meshes of Pig 2. In this case, the registration improves widely the correspondence between the post-and the pre-operative mesh as suggested by the point-to-mesh errors (Table 7.1).

Mesh from Pre-op image Before registration

After registration Figure 7.3: Comparison between the pre-operative arterial mesh in blue and the post-operative meshes (in red: before the registration and in green: after the registration) of Pig 2. The correspondence between the post-and the pre-operative meshes is improved.

Evaluation of the Smooth Vessels Trees

In order to evaluate the creation of the smooth vessels trees, DICE, Positive Predictive Value (PPV) and sensitivity scores between the original and the smoothed vessel images are computed. The values for Pig 4 are reported in Table 7.2. They show a good correlation between the original segmented vessels and the generated ones. Given V m and V s , the volume of measured (respectively simulated) necrotic area, then the DICE score (DICE), the sensitivity (S) and the Positive Predictive Value (P P V ) are defined as:

DICE = 2|V m ∩ V s | |V m | + |V s | S = |V m ∩ V s | |V m | P P V = |V m ∩ V s | |V s | Figure 7
.4 illustrates the inconsistent areas, which explained the errors in DICE, PPV and sensitivity for Pig 4 (Table 7.2).

Portal vein after smoothing

Portal vein before smoothing

Arteries after smoothing

Arteries before smoothing Figure 7.4: (Left): Comparison between the arterial meshes before and after the smoothing process of Pig 4. Some vessels, which have a strong stenosis are canceled after the smoothing process. (Right): Comparison between the portal vein meshes before and after the smoothing process. Irregularities in the segmentation are smoothed as illustrated by the green circles.

Verification of the Parameter Estimation Framework

Synthetic Data Generation

In order to verify the parameter estimation framework, we consider a synthetic case on a regular cuboid domain to speed up the process (Figure 7.5). As we try to estimate the heat conductivity and the heat capacity, only the diffusion is considered in this case.

Personalized parameters found in Section 7.8 are used. We apply the typical clinical RFA protocol: during three minutes, we increase linearly the temperature at the tips of the probe, then we heat at 105 • C for six minutes, and finally no temperature is imposed for 3 more minutes as illustrated on Figure 7.6, right panel. With this protocol, a necrotic area, the delivered power and the temperature during the cooling are simulated as shown on Figure 7.6. Those data are then used as ground truth. The main parameters: dt and c t are estimated by minimizing the cost function defined in Eq. 6.1. The range of parameters values used [Altrogge 2012] are reported in Table 7.3 and the other parameters are fixed to the nominal values, reported in Table 6.1. The temperature distribution over ablation on the synthetic case. During the heating period, the temperatures are imposed at the nine tips of the probe and the delivered power is simulated. During the cooling period, the temperatures are simulated.

X-axis

Evaluation

We compare two gradient-free optimization methods available in DAKOTA1 : the Constrained Optimization BY Linear Approximations (COBYLA) and the pattern search method (PS). Those methods are sequential trust-region algorithms. Initially, the total domain of parameters is visited (Table 7.3) and then the region is contracted. Using COBYLA, we manage to obtain the estimated parameters with 0.8% error on dt and 0.6% on c t in 30 minutes after 22 iterations with an error of 7.8 10 -5 of the cost function and a mean of the symmetric point-to-mesh error of 10 -4 mm. Similarly, using PS, we manage to obtain the estimated parameters with 0.4% of error on dt , 0.3% on c t in 154 minutes after 117 iterations with an error of 7.2

Notation Parameter Name Min -Max dt heat conductivity 0.256 0.768 W(m K) -1 c U t heat capacity of undamaged tissue 18 54 J(kg K) -1 Table 7.3: Ranges of parameters values explored in the optimization framework.

10 -6 of the cost function and a similar mean of the symmetric point-to-mesh error of 10 -4 mm. The convergence curves of both algorithms are shown on Figure 7.7.

We choose the COBYLA method to perform the personalization on the pigs as it requires less iterations of the forward model to estimate the parameters with similar accuracy as PS.

Sensitivity Analysis of the Model Parameters

In order to reduce the number of parameters to estimate for the personalization of the model, a sensitivity analysis on a synthetic case is performed. We use the same synthetic case as the one used in the previous section 7.3 and described in Figure 7.5. During the heating phase, the temperature at the tips of the probe is imposed as it is shown in Figure 7.6. For this reason, we chose to show only the cooling phase of the temperature. To assess the effect of one given parameter, we manually change its value and run different simulations, all other things remaining equal. Figure 7.8 illustrates the effect of the heat capacity of undamaged cells, c U t on the delivered power and on the cooling temperature distribution, especially on the time of the peak value of the delivered power. Figure 7.9 illustrates the effect of the heat capacity of vulnerable cells c V t on the delivered power and on the cooling temperature distribution. This parameter as an effect both on the delivered power and on the cooling temperature distribution. 7.10 illustrates the effect of the heat capacity of necrotic cells, c N t on the delivered power and on the cooling temperature distribution. This parameter does not largely affect the delivered power neither the cooling temperature.

Finally, Figure 7.11 illustrates the effect of the heat conductivity, dt on the delivered power and on the cooling temperature distribution. The heat conductivity has a clear impact on the cooling temperature distribution. We chose to estimate the heat conductivity, dt for its effect on the cooling temperature and the heat capacity of undamaged cells c U t for its effect on the time of the peak value of the delivered power.

Evaluation on Synthetic Data

We want to verify the use of our parameter estimation framework in the case where the considered liver is not healthy. To do so, we use a synthetic case on a regular cuboid domain as the available data come from healthy pigs. We consider the case of a cirrhotic liver, where it is known that there is an "oven effect". The propagation of the temperature is more difficult in a cirrhotic liver than in a healthy liver, the cirrhotic tissue behave like a thermal insulator, preventing heating outside the tumor targeted for the ablation [START_REF] Livraghi | [END_REF]]. In order to mimic this effect, the heat conductivity is divided by two. With this value, we apply the typical clinical RFA protocol, and we obtain a "cirrhotic" necrotic area, the delivered power and the temperature during the cooling phase, which differs from those obtained with the nominal value, as illustrated on Figure 7.12, right panel. Those data are then used as ground truth in order to estimate the main parameters: dt and c t by minimizing the cost function defined in Eq. 6.1.

Simulated lesion with nominal parameters

Simulated lesion with "cirrhosis" parameters

Figure 7.12: (Left): The parameter estimation framework proposed is able to perfectly find the "cirrhotic" heat conductivity. (Right): The lesion obtained with the "cirrhotic" heat conductivity in green is enclosed within the lesion obtained with the nominal conductivity coefficient in blue.

The parameter estimation framework is used, and we are able to find exactly the cirrhotic value of the heat conductivity as illustrated on Figure 7.12, left panel. We show on a synthetic case, that our method allows to find "cirrhotic " value of the heat conductivity in the case of a diseased liver. We can now evaluate our method on the available swine data, whose liver tissue are healthy.

Evaluation on Swine Data

Systematic Study

During our experiments, the measured electrical power appears consistently to reach its maximum during the increase in temperature, before reaching the plateau of 105 • C (Figure 7.13, left panel) in all cases. and the four remaining temperatures are imposed from a linear interpolation from those measurements (four green curves). After the use of personalized parameters, the error between the measurement and the computation is reduced from 18.9 to 8.8 W for the power and from 3.5 to 3.0 • C for the temperature. A zoom on the cooling phase is presented in the following Figure 7.14.

This phenomena cannot be explained by a constant heat capacity which would lead to a peak after the plateau is reached. Instead, this observation suggests that the cells reach their vulnerable state faster. Thus we update our necrosis model accordingly: after studying the ODEs of the model and a sensitivity analysis on its parameters non reported here, δ is modified such that tissues reach very fast their vulnerable state which entails a significant change of heat capacity (Table 7.4).

The model is evaluated on twelve ablations performed in five swines (several surrogate tumors ablated for some pigs). To understand the effect of the perfusion on the necrotic lesion, the simulations are performed first without the advection term. We only take into account the cooling effect of the blood vessels, and we will account for the hepatic perfusion later by adding the advection term in a second step (Section 7.7.1).

The applied RFA protocol is not exactly the same for all ablations. Eight ablations are performed through several short cooling and heating periods, whereas the other four ablations included only one long final cooling stage after a continuous cooling Figure 7.14: Zoom on the cooling phase of the temperature distribution after ablation from tumor 4-3. The personalized parameters do not affect the cooling temperatures.

heating period (it is the case for tumor 4-3 as shown on Figure 7.13). For all pigs, nominal values of parameters (reported on Table 6.1) are employed. In each case, the simulated lesion was compared with the registered ground-truth. Quantitatively, symmetric point-to-mesh errors implemented like in [START_REF] Zheng | [END_REF]] computed between the simulated lesion and the registered post-operative groundtruth showed good prediction of the necrotic extent (average over all ablations of 5.32 ± 3.62 mm of mean point-to-mesh error) as illustrated in red on Figure 7.16. Similarity scores are also computed and shown on Figure 7.17 in red: average over all ablations: DICE of 44.2%, Sensitivity of 47.0%, and PPV of 52.7%. We achieved errors of 25.9 W and 5.1 • C on average between measured and simulated values of cooling temperature and delivered power (Figure 7.18 in red). Values computed with nominal parameters (nominal) and without the advection term are reported in Table 7.10 for Pig 1 and Pig 4, Table 7.11 for Pig 2, Table 7.12 for Pig 3 and Table 7.13 for Pig 5. The prediction of the necrosis extent was valid up to 5mm on average which can be considered as sufficient for clinical applications. Qualitatively, as one can see on Figure 7.19 (red curves), the simulated heat power and temperature were close to the heat power and the temperature given by the RFA probe itself. During the heating phase, the temperature is imposed, the power is simulated and compared to the measured one. On average, the error is reduced after the use of personalized parameters: from 35.1 W to 29.8 W. It is mainly reduced on the first heating period (from 50.3 W to 26.5 W). During the cooling phase, the simulated (non-imposed) temperature is compared to the measured one, whenever measurement is available (the RFA power was turned off during the cooling phases). On average, the error is also reduced: from 5.9 • C to 1.9 • C. During the heating phase, the temperature is imposed, the power is simulated and compared to the measured one. The error in term of delivered power is constant (35.0 W versus 35.1 W). During the cooling phase, the simulated (non-imposed) temperature is compared to the measured one, whenever measurement is available (the RFA power was turned off during the cooling phases). There is no real difference in term of cooling temperature error: 6.3 • C versus 5.9 • C.

DICE SENSITIVITY PPV

Towards Personalization

As illustrated on the pipeline of Figure 6.2, the framework enables parameter estimation too. COBYLA is used to minimize the cost function (Eq. 6.1) as only a few forward simulations (typically 20) are required as described in Section 7.3. The data came from healthy pigs of similar age and weight, we hypothesize that the parameters would be the same for all of them. The vulnerable rate coefficient is also adjusted to match the raise in delivered power as detailed in Section 7.6.1. As the effect of the hepatic perfusion has not been clearly identified so far, we do not consider the advection in the parenchyma in this part of the study. We run the simulation without taking into account the effect of the blood flow in the parenchyma, but only the heat sink effect of the blood flow in the large visible vessels. The results of the personalization are then compared with simulations where the advection term is not accounted as well: (perso) versus (nominal). This is a proof of concept not only for the ability to personalize the model but also for its predictive power by evaluating the simulation results on ten different tumors. Briefly, the minimization of the error between measured and simulated values of power and temperature was done only on two tumors with a long final cooling stage (2 different tumors: Pig 1 and Tumor 4-3 of Pig 4) as it was long enough to observe reliably the effect of the conductivity dt , yielding two sets of personalized values. In both cases, the values independently found were really close: the same value of heat capacity c U t was estimated, and the conductivity values were almost equal to the nominal value, as expected. The values are reported on Table 7.4 and shown on Figure 7.21. The estimated heat capacity and the nominal conductivity were then used to simulate RFA on the ten remaining cases and errors in temperature, heat power and necrosis size were evaluated (Table 7.10 for Pig 1 and Pig 4, Table 7.11 for Pig 2, Table 7.12 for Pig 3, Table 7.13 for Pig 5, (perso)).

Small errors were obtained in those cases too, without previously having fit the parameters for those tumors. With those changes, the simulated electrical power matches much better the measured one in terms of time to peak: average error between measured power and simulated power reduces from 25.9 W to 20.4 W on average (Figure 7.18 in green, left panel). However the mean point-to-mesh error (5.99 mm on average) and the error in temperature (7.8 • C on average) do not change significantly (in green in Figure 7.16 and Figure 7.18,right panel) as well as the similarity scores: Dice of 42.2%, Sensitivity of 43.0% and PPV of 54.6% on average (Figure 7.17 in green). Figure 7.22 and Figure 7.13 illustrate the improvement in term of delivered power after personalization on Pig 1 and tumor 4-3 (red curves versus green curves).

Although the personalization has only been performed on two pigs, we can see this improvement on tumor 3-1. Figure 7.19 shows how the first cooling phase is closer to the measurement given by the RFA probe itself with personalized parameters (green curve) rather than with the nominal parameters (red curve). The heat capacities of the vulnerable neither the necrotic cells (c V t , c N t ) were modified, thus the simulations with personalized parameters and with nominal parameters of the delivered power for the following heating phases are similar. As the estimated value of dt after personalization is the nominal value, the personalized parameters do not have any impact on the cooling temperature, not on the tumors where the personalization has been performed (Figure 7.23), neither on the others (Figure 7.19, cooling heating Figure 7.22: Personalization of the heat capacity and the conductivity using (Left) the delivered power curve over ablation (Right) the temperature distribution after ablation from tumor 1. The error between the measurement and the computation is reduced from 19.9 to 7.5 W for the power, whereas it is constant for the cooling temperature (11.1 • C versus 10.3 • C). A zoom on the cooling phase is presented on the following Figure 7.23. Down). We do not have the measured cooling temperatures for this tumor 3-1, however we can see that this value of dt give accurate simulations, since there is no jump between the computed temperatures at the end of a given cooling period and the measured temperatures at the beginning of the following heating period for the different heating/cooling stages at times: t = 300s, t = 480s, t = 680s and t = 890s (Figure 7.19, Down). In order to quantify the effect of hepatic perfusion, the simulations were performed by adding the advection term for all pigs whenever possible. In some cases (Pig 2 and Pig 4), the segmentation of the inlet of vena cava was not possible, thus the impossibility to get results from the CFD solver. Nominal values of parameters (reported on Table 6.1) were employed as well. The advection has an effect on the shape of the lesion, however in this case, it does not improve the point-to-mesh error: 6.9 ± 5.2 mm versus 4.8 ± 4.4 mm.

cooling heating

Sensitivity Analysis on Swine Data

Hepatic arteries

The quantitative results are reported in Table 7.10 for Pig 1, Table 7.12 for Pig 3 and Table 7.13 for Pig 5. (nominal) stands for computation without the advection term and (advection) with the advection term. On average, the effect of the advection cannot be clearly highlighted. The mean of the point-to-mesh error over the 6 ablations computed with advection is 4.96 ± 3.50 mm (Figure 7.16 in blue), the average DICE score is 45.6%, the average sensitivity is 41.8% and the average PPV is 53.4% (Figure 7.17 in blue). On average, the error on the simulated power is 22.1 W and 4.7 • C on the temperature (Figure 7.18 in blue). Those values are comparable to the values obtained without the advection term. However, even if it is not the case globally, the advection has an effect on the computed lesion as well as the computed temperature and power depending on the ablation location. In some cases, the advection reduces the error in term of necrotic extent and temperature distribution (tumor 5-1), however it is not true in every case (Pig 1), but the advection term do change the shape of the necrotic extent (Figure 7.24). Quantitatively, this effect of the advection on the lesion shape can be seen on Table 7.14 for Pig 1 and Pig 3, Table 7.15 for Pig 5. In the cases where the advection has an effect, the sphericity of the lesion decreases and the maximum diameter of the lesion increases. The advection has an effect on the temperature distribution at the tips of the probe as well, depending on the ablation location (Figure 7.25). cooling heating Figure 7.27: Zoom on the cooling temperature distribution after ablation for tumor 1. When the advection term is taken into account, the cooling temperatures at the different tips of the probe do not have to same cooling speed (blue curves). Whereas, they cool down similarly without the advection term (red curves). run without advection. In this case, it could be explained by the fact than the perfusion tends to cool down the parenchyma and thus more power is needed to maintain the targeted temperature. In the case of short heating and cooling stages, we do not observe this phenomena (Figure 7.20), it might be because the liver temperature drops quickly when the delivered power is decreased.

Tumor Probe

Post-op lesion

Effect of the Registration

One part of the evaluation of the computation is done by a comparison between the simulated necrotic extent and a lesion segmented on a post-operative image registered to the pre-operative image. The modality of the post-operative image is not always the same (Table . 7.5). Table 7.5: Modalities of the pre-and post-operative images used. The segmentations of liver, vessels, tumor on the pre-operative image and of necrosis on post-operative image do not come from the same modality for each pig.

Pig

pre-op image post-op image

1 CT CT 2-1 CT T1 + gad 2-2 CT T2 3-1 CT-Venous N.A 3-2 " T2 3-3 " " 4-1 CT CT 4-2 " " 4-3 " " 5-1 CT-Port T2 5-2 " " 5-3 " "
In some case, the resolution and the size of the post-operative image do not allow an accurate registration of the lesion. It is the case for tumor 2-2 for example (resolution: 0.78x0.78x4.8 mm, size: 384x384x36 mm). To check the sensitivity of the computation to the registration, a new registration of the post-operative lesion is performed for this tumor. We assume that the barycenter of the post-operative lesion is the same as the barycenter of the simulated necrotic extent, and we rigidly translate the post-operative lesion to align the two barycenters (Figure 7.28,right panel).

With this simple registration method, the results are significantly improved (Table. 7.6). The lesion registration strongly affects the error measured between the computed necrotic area and the ground truth. 5.03, 0.11, 23.30, 8.37, 16.59 13.8 9.8 23.3 (nominal-0.1) 8.03, 4.85, 0.10, 22.28, 7.83, 15.58 16.8 13.1 23.5 (perso) 8.44, 5.00, 0.10, 22.99, 8.27, 16.35 14.3 10.4 22.6 translation (nominal-0.5) 3.87, 2.95, 0.10, 18.00, 2.98, 10.78 57.3 40.7 96.5 (nominal-0.1) 3.12, 2.65, 0.11, 17.08, 2.31, 9.38 66.3 51.7 92.1 (perso) 3.63, 2.87, 0.14, 17.68, 2.70, 10.24 60.1 43.9 94.9

Tumor

Simulated lesion

Post-op lesion after non-rigid registration

Post-op lesion after translation

Effect of the Probe Position

The position of the probe is known from an intra-operative anatomical MRI. In some cases, the resolution and the image size do not allow an accurate registration of the probe. It is the case for tumor 3-2 for example (resolution: 1x1x3.5 mm, size: 320x250x30 mm).

A new computation is performed for this tumor in order to check the sensitivity of the computations to the probe position. The same configuration was used except that we manually put the probe inside the registered post-operative necrotic area (position 2 in Figure 7.29). With this probe position, the results are significantly improved (Table. 7.7). The probe position affects also the error measured between the computed necrotic area and the ground truth (Figure 7.29). It is equivalent to a registration error since by moving the probe, we move the computed necrosis.

Cell Death Model Evaluation

As the parameters of the cell death model have been modified after parameter estimation, the accuracy of the model has to be verified. We compare the simulated necrotic extent for tumor 2-2 with the results of a lethal thermal dose model, often used in clinical settings [Sapareto 1984] (Figure 7.30). The thermal dose at each voxel is calculated according to this empirical model:

t 43 = tc t=0 C (43-T ) ∆t (7.1)
where t 43 is the thermal dose in equivalent minutes at 43 • C, T is the average temperature during the time interval ∆t, and C is a constant resulting from the Arrhenius model.

C = 0.5 when T ≥ 43 • C C = 0.25 when T < 43 • C
A thermal dose of 240 (in equivalent minutes at 43 • C) is considered lethal. Quantitatively, both models give similar point-to-mesh errors and DICE scores (Table . 7.8), however the volume and maximum diameter of the necrotic extent are larger with the thermal dose model (Table. 7.9). The thermal dose model gives really spherical necrotic extents (Table . 7.9, sphericity values closer to 1 for thermal dose model). A more throughout study is necessary to understand the differences between those two models of cellular necrosis. However, the proposed three-state model enables a spatial and time-varying heat capacity (c U t ,c V t ,c N t , for undamaged, vulnerable and necrotic cells). We saw in Section 7.4 that those parameters affect the delivered power and the cooling temperature. With the thermal dose model, only two states of the cells could be considered: they are either undamaged (if T D < 240) or necrotic (if T D ≥ 240). Thus only two values of the heat capacity could be used. Thermal Dose (nominal-0.5) 8. 25, 5.71, 0.03, 21.1, 7.36, 18.87 13.9 23.2 9.9 (nominal-0.1) 9.05, 6.50, 0.06, 24.99, 7.76, 21.10 14.3 29.9 9.4 (perso) 8.68, 6.13, 0.15, 22.77, 7.61, 20.09 13.9 25.9 9.5 Cell Death (nominal-0.5) 8. 52, 5.03, 0.11, 23.30, 8.37, 16.59 13.8 9.8 23.3 (nominal-0.1) 8.03, 4.85, 0.10, 22.28, 7.83, 15.58 16.8 13.1 23.5 (perso) 8.44, 5.00, 0.10, 22.99, 8.27, 16.35 14.3 10.4 22.6

Discussion

We have presented a multi-physics model based on multi-modal medical images, evaluated on pre-clinical data. We rely on LBM to solve the bioheat equations and to compute the blood flow in the liver. Despite possible biases in the probe location and in establishing correspondences between the post to the pre-operative images due to registration errors, our model provided new insights, and we show that parameter estimation is possible to reduce the bias introduced by the use of nominal parameters. This opens new opportunities for RFA planning and guidance.

Model Limitations

The impact of the segmentation process has to be taken into account as well as the model limitations (assumptions: we assume that the measured power is actually the heat power through radio-frequency, with no loss, etc).

We are not able to accurately locate the position where we measure the pressure, but we assume that there is no variation in pressure in a small neighborhood of the two outlets.

The Neumann Boundary condition applied on the border of the liver can be discussed. If the ablation location is subcaspsular (tumor 5-1, tumor 4-1, tumor 3-3), the volume of the necrotic area is larger than in the other cases (Table 7.14, Table 7.15), as the heat does not dissipate through the liver capsule but reflects on the border and increases artificially the heat diffusion inside the liver as it is illustrated in Figure 7.31. However, the liver has a high regulation capacity thanks to its rich vascular system. If the heating is applied close to the border of the liver, it is not clear how the heat will dissipate, whether or not it will affect the neighboring organs.

Finally, during the RFA computation, the temperature at the tips of the probe is imposed. It will be more accurate to directly impose the delivered power and thus simulate the cooling as well as the heating temperature distribution. It should be done in the future, but going from a Dirichlet boundary condition on the temperature at the tips of the probe to a Dirichlet boundary condition on the delivered power is not straightforward.

Effect of the Registration

It should be noted that the registration method used in this paper differs from the registration used in [Audigier 2015b]. In [Audigier 2015b], due to the uncertainty in the registration of the post-operative image to the pre-operative image, the necrotic lesion segmented on the post-operative image was registered rigidly to the preoperative image by aligning its barycenter with the barycenter of the simulated necrosis which could introduce bias in the analysis. The impact of those registration processes has been evaluated on one tumor in this chapter.

Towards Personalization

In clinical RFA of liver tumors, we cannot assume that the biophysical parameters are the same for all patients, as assumed here for the five pigs, and parameter personalization is therefore required. Here, we propose to adjust model parameter for a better evaluation of RFA, based on the delivered electrical power during ablation and the temperature drop during cooling. Since our approach does not rely on ablated regions, the parameter estimation leads to an increased confidence in the computed temperature, in the case where no temperature map is available. By fitting the temperature and the delivered heat power, the novel approach can estimate the temperature around each tumor at any time during the ablation. This additional information could be used as surrogate to assess the amount and location of damaged tissue during the intervention (cells receiving excessive heat but without being necrosed) surrounding the ablated region.

Furthermore, the probe temperature and delivered power are information that are readily available, in real-time, from the RF system and therefore could be used for therapy guidance. Even if the computation time does not allow an interventional set-up yet, it is a first step towards the personalization of RFA computation.

The extent of the necrotic area is mainly controlled by the conductivity dt and the heat capacity of vulnerable cells c V t (not of undamaged ones c U t ). Since dt was not changed after optimization, there was no significant difference in terms of necrotic area after personalization (mean of the point-to-mesh error of 5.99 mm versus 5.32 mm on average) or cooling temperature (error of 7.8 • C versus 5.1 • C), despite a better match for the measured power (error of 20.4 W versus 25.9 W) as illustrated on Figure 7.19 and Figure 7.22, and reported in Table 7.10, Table 7.11, Table 7.12, Table 7.13. Moreover the point-to-mesh errors were of the order of 5mm; it suggests that simulations with optimized dt and c U t are realistic in terms of necrotic area, power and temperature predictions which was the objective. These results confirmed the stability of the personalization framework.

Conclusion

In this chapter, we evaluated our model on a pre-clinical data. The approach is evaluated on five swines and twelve ablations. It gives promising results, but more importantly provides several lessons learned.

Despite a very comprehensive model and a complete pre-clinical study for its validation, many errors are involved, at the modeling level, but also especially at each step of the pre-processing. We manage to identify most of those errors thanks to the pre-clinical study. It would not have been possible on clinical data. The use of invasive catheter to measure the pressure could not be done on a clinical study, neither the acquisition of Phase-Contrast MRI, intra-operative MRI with the probe, injected CT scans with the three different phases, which are not part of the standard clinical workflow for RFA.

Even if model personalization is premature since the computation time has to be reduced first, we show that it can be done based on intra-operatively simple quantities (the probe temperature and delivered power) since we rely on a LBM implementation on graphics processing units (GPU): a forward simulation is faster than real time. As a first step towards personalization, we evaluate the discrepancy in terms of temperature and delivered power, and we show that key biophysical parameters can be estimated leading to promising predictions. In clinical settings, due to the large variety of diseases treated by RFA (cirrhosis, fibrosis, etc), the proposed method should be suitable to get patient-specific parameters from easily accessible data.

A study of the errors is a necessary next step before using model-based therapy in the clinical routine. A study of how the pre-processing errors (due to segmentation, post-to-pre registration, probe registration, vessels smoothing) propagate to the final computational outcome will be of great interest. It will help to understand the modeling errors as well and be more confident on the resulting necrotic lesion for example. The evaluation could be improved by the acquisition of images of better resolution for the intra-operative MRI with the probe for example, but taking into account that with the current framework, the acquisition time is limited by the well-being of the pig. The acquisition of the Phase-Contrast MRI could be improved as well. Moreover, even if today MR thermometry does not provide such a good resolution (Figure 7.32 shows MR thermometric images acquired during the ablation of tumor 4-1), in the future, it could add more information on the spatial extent of the temperature to further improve the validation.

-60 °C 60 -75 °C 75 -88 °C

Before Ablation

During Ablation To resume, this thesis presents novel contributions in RFA computational simulation, which can be summarized as follows:

1. A complete patient-specific geometry including hepatic venous and arterial circulation system;

2. The patient-specific modeling including simultaneously the cooling effect of large vessels and of the perfusion within the parenchyma;

3. The discretization method (LBM), which is fast and verified against an analytical solution.

Model Evaluation on Clinical Data

A necessary step before deploying model-based therapy planning or guidance in clinical settings is a clinical validation with extensive data on larger populations to evaluate the computational model of RFA. We showed that it will be possible due to the level of integration of the proposed framework. As a first step towards clinical application of a RFA computational model, this thesis presents a clinical evaluation on a database of ten patients of the developed RFA computational simulation.

The results we obtained by comparing the real necrotic lesions from post-operative images with the simulated lesions suggest a promising correlation. The effect of the blood perfusion inside the liver parenchyma as well as the blood flow through the large vessels has also been highlighted. Finally, a first experiment of tissue parameter adjustment on one patient indicates the need of model personalization.

Our system constitutes a first step towards the clinical validation of RFA modeling, which can be summarized as follows:

1. The evaluation of the proposed RFA model on a dataset of ten patients;

2. Experiments testing several hypotheses in the discussion: effect of the probe position, effect of the segmentation;

3. The effect of the parenchyma and venous blood flow on the resulting necrotic lesion (ablated tissue);

4. The importance of the personalization of the model biophysical parameters.

Towards Model Personalization

In clinical RFA of liver tumors, due to the large variety of diseases treated (cirrhosis, fibrosis, etc) we cannot assume that the biophysical parameters are the same for all patients, and parameter personalization is therefore required. We showed that the proposed method could be suitable to get patient-specific parameters from easily accessible data. As we rely on the Fully Coupled LBM CFD solver, no advanced meshing techniques are required. All the computations are directly done from patient images using LBM: heat propagation and cell death modeling as well as the heat sink effect of blood vessels and porous circulation in the liver. Since we rely on an implementation on graphics processing units (GPU): a forward simulation is faster than real time. Two personalization strategies were tested and evaluated on clinical and pre-clinical data.

First, an automatic estimation of the main parameters of the model is proposed and evaluated on clinical data. In this case, the optimization algorithm minimizes the error between the ground truth and the computed necrotic lesion. Therefore, the method needs several tumors for validation. Promising results are achieved with the use of patient-specific parameters on three patients and seven tumors ablated.

Then, a second personalization strategy was presented. In this case, the optimization algorithm minimizes the error between the computed values of the probe temperature and delivered power and the actual values measured by the RFA device. It is evaluated on a pre-clinical dataset of five swines. The personalization of the sensitive tissue parameters has been performed on two surrogate tumors and its predictive power evaluated on the remaining ten surrogate tumors. We showed that key biophysical parameters can be estimated leading to promising predictions. This method is based on intra-operative simple quantities (the probe temperature and delivered power) that are readily available, in real-time, from the RF system and therefore could be used for therapy guidance.

Even if the computation time does not allow an interventional set-up yet, we showed that model personalization can be done based on intra-operative quantities. It is a first step towards the personalization of RFA computation.

Pre-clinical Model Evaluation

The proposed model has been evaluated on pre-clinical data from a comprehensive experimental set-up specially designed for radiofrequency ablation (RFA) model validation. To the best of our knowledge, it constitutes the largest and most comprehensible in vivo study reported so far. The model considered is complete and takes into account the main biophysical mechanisms (heat transfer, cellular necrosis, hepatic blood flow) involved in RFA. The pre-clinical dataset includes five swines and twelve ablations, leading to promising evaluation results, but more importantly to several lessons learned. Our pre-clinical study was challenging to establish due to the complexity of the experimental protocol. Many images from different modalities and at different times of acquisition were involved. Before accounting for all the information from the images into the RFA model, an impressive amount of preprocessing was required: segmentation, vessel smoothing, probe modeling, probe registration, image registration. Moreover, the setting of subject-specific boundary conditions for the blood flows and blood pressure was not straightforward. Many errors are then involved, coming from the model as well as the pre-processing. However, thanks to the pre-clinical dataset used, we manage to identify most of them, that will not be possible with a clinical study. The effects of the segmentation, the probe registration, the post-operative lesion registration and the vessel smoothing have been evaluated. At the modeling level, the effects of the advection, the probe position, the cell death model as well as the sensitivity of the model biophysical parameters were highlighted.

Perspectives

Model Improvements

Although the implemented model includes the main important biophysical mechanisms involved in radiofrequency ablation of hepatic tumors (heat transfer, cellular necrosis, hepatic blood flow), several improvements can be performed. For instance, a weak coupling between the blood flow and the heat transfer model is currently considered. The blood flow has an influence on the temperature distribution through the WK model and through the reaction term in the Pennes model but the temperature does not affect the blood flow. This weakly coupling model allows us to speed up the calculations since the blood flow distribution is computed only once, at the beginning of the simulation. But a full coupling of blood flow and heat transfer models should give more accurate predictions, and enable us to consider coagulation and water evaporation, thus improving on the presented model.

Furthermore, we realized through the pre-clinical study that the delivered power is information that can be available easily in real-time from the RFA probe itself. In the current implementation of the model, the probe temperature is imposed through Dirichlet boundary conditions during the heating phase. We are then able to simulate completely the delivered power but to simulate the probe temperature during the cooling phases only. By imposing the delivered power instead of the probe temperature, the model will be improved and we will be able to simulate the probe temperature during the total duration of the procedure. Imposing Dirichlet boundary condition on the power in the current LBM implementation is not straightforward, but it is worth investigating this aspect. It would be a more coherent way of simulating RFA.

Finally, some interesting aspects to improve the ongoing RFA intervention have not been considered in this work, but the current RFA model opens possibilities to an automatic estimation of the optimal locations of the probe and the electrical power to deliver in order to achieve a proper tumor ablation for example.

Model Validation

A study of the modeling errors is a necessary next step before using model-based therapy in the clinical routine. Understanding how the pre-processing errors (due to segmentation, post-to-pre registration, probe registration, vessels smoothing) prop-agate to the final computational outcome is necessary. It will help to make the model more realistic and to have a confidence measure on the resulting necrotic lesion for example. Especially, we saw that the segmentation has a high impact on the computed lesion as the blood flow computation through the CFD solver is strongly correlated with the vessel segmentation. The post-to-pre and the probe registrations influence the validation of the model output. Since modeling of RFA is still an open challenge, the present RFA model is a valuable asset, and future directions include a more throughout validation on a larger cohort of patients.

Moreover, the validation could be improved by the acquisition of good resolution and reliable MR thermometric images. As data acquisition techniques improve quickly, we believe that it would be possible in the near future. It could add more information on the spatial extent of the temperature even far from the RFA probe and thus further improve the validation.

Model Personalization

Model personalization is premature since the computation time has to be reduced first. But we showed on clinical data that it can be done in the case of multiple tumors inside the liver. Adaptation for RFA under image-guidance using intraoperative image and no post-operative image has been considered. However, a necessary step before deploying this method in clinical settings is a pre-clinical validation to evaluate the computational model of RFA and also to consider potential safety issues of the proposed application.

By evaluating the discrepancy in terms of temperature and delivered power on pre-clinical data, we showed that intra-operative simple quantities (the probe temperature and delivered power) can be used for personalization of model-based therapy. Moreover, since this approach does not rely on ablated regions, the parameter estimation leads to an increased confidence in the computed temperature, in the case where no temperature map is available. It can estimate the temperatures around each tumor at any time during the ablation. This additional information could be used as surrogate to assess the amount and location of damaged tissue during the intervention (cells receiving excessive heat but without being necrosed) surrounding the ablated region, helping the clinical monitoring during the RFA intervention.

Conclusion

Several research groups have already proposed different frameworks for the modeling of radiofrequency ablation (RFA). They have showed the potential benefits of such models but they have also pointed out the challenges involved. Through this thesis study, we showed that important advances have been made in the patient-specific modeling of RFA of hepatic tumor. The results we obtained are promising and after further future improvements, we believe that such personalized model could be used in clinical environment. Although our research was focused mainly on the RFA of hepatic tumor, the same pipeline of methods can be used for tumors located in other organs, and could be adapted easily to other percutaneous ablative techniques like High Intensity Focus Ultrasound, Microwave Ablation or Laser Ablation for example, according to the properties of the current ablation procedure considered. However, more research is required to tackle the issues we faced at each step of the work flow (segmentation, registration, model validation and so on...). The study and validation of model simplification will be necessary before its integration in the medical practice. As data acquisition techniques and our understanding of the biophysical mechanisms involved in RFA improve quickly, such models will become more accurate and more realistic in the near future. We believe that patient-specific simulations will become a common tool in the therapy planning and guidance for a specific patient.

Parameter Estimation For Personalization of Liver Tumor Radiofrequency Ablation.

In MICCAI Workshop on Abdominal Imaging-Computational and Clinical Applications 2014 (pp. 3-12). A statistical description of a system can be explained by the distribution function (probability of finding particles within a certain range of velocities at a certain range of locations at a given time) f (r, c, t) where f (r, c, t) is the number of molecules at time t positioned between r and r + dr which have velocities between c and c+dc [Mohamad 2011]. An external force F will change the velocity of the molecule from c to c + Fdt and its position from r to r + cdt. The number of molecules f (r, c, t), before applying the external force is equal to the number of molecules after the disturbance f (r + cdt, c + Fdt, t + dt) if no collisions take place between the molecules:

f (r + cdt, c + Fdt, t + dt)drdc -f (r, c, t)drdc = 0
However, collisions can take place between the molecules. The rate of change between final and initial status of the distribution function is called collision operator, Ω. Hence:

f (r + cdt, c + Fdt, t + dt)drdc -f (r, c, t)drdc = Ω(f )drdcdt
Dividing the above equation by drdcdt and as the limit dt → 0, yields where a is dc dt , the acceleration and can be related to force F by Newton's law:

df dt = Ω(f )
a = F m
Therefore, the Boltzmann transport equation can be written as:

∂f ∂t + ∂f ∂r c + ∂f ∂c F m = Ω(f )
The relations between the above equation and macroscopic quantities such as fluid temperature T , fluid velocity vector u are:

T (r, t) = f (r, c, t)dc T (r, t)u(r, t) = cf (r, c, t)dc

A.1.2 The BGKW Approximation

It is difficult to solve Boltzmann equation because the collision term is very complicated. Bhatnagar, Gross and Krook (BGK) in 1954 introduced a simplified model for collision operator:

Ω(f ) = 1 τ (f eq -f )
where τ is called the relaxation factor. The local equilibrium distribution function is denoted by f eq , which is Maxwell-Boltzmann distribution function.

After introducing BGKW approximation, the Boltzmann equation (without external forces) can be approximated as:

∂f ∂t + c.∇f = 1 τ (f eq -f )
In Lattice Boltzmann Method, the above equation is discretized and assumed valid along specific directions. Hence, along a specific direction, the discrete Boltzmann equation can be written:

∂f i ∂t + c.∇f i = 1 τ (f eq i -f i )
This equation can be discretized as:

f i (r + c i ∆t, t + ∆t) = f i (r, t) + ∆t τ [f eq i (r, t) -f i (r, t)]
Let us consider: w = ∆t τ Then:

f i (r + c i ∆t, t + ∆t) = f i (r, t) + w[f eq i (r, t) -f i (r, t)]

A.1.3 Equilibrium Distribution Function

The key element in applying LBM for different problems is the equilibrium distribution function f eq . For particles moving in a medium with macroscopic velocity u, the normalized Maxwell's distribultion function can be written as:

f eq = T 2π/3 e -3 2 (c-u) 2
which can be written as,

f eq = T 2π/3 e -3 2 c 2 e 3(c.u-u 2 /2)
where c 2 = c.c and u 2 = u.u.

Using Taylor series expansion for e x and expanded around the stationary state:

f eq = T 2π/3 e -3 2 c 2 [1 + 3c.u - 3 2 u 2 + ...]
Let us write the general form of the equilibrium function as:

f eq i = T ω i [A + Bc i .u + C(c i .u) 2 + Du 2 + ...] where A, B, C, and D are constants and need to be determined. T stands for scalar parameter (in our case: the fluid temperature), which is equal to summation of all the distribution functions, i.e:

T = N i=1 f eq i
where N is the number of lattice links.

In our case, we will just consider: The last term can be written as: w i e i,α e 2 i,β u α = 0

f eq i = T ω i [A + Bc i .
The third momentum is:

T c 2 c 2 s = N i=0 f eq i c 2 e i 2
This could also be written in a matrix, since: ( t e i .e i ) α,β = M k=0

( t e i ) α,k (e i ) k,β = e i,α e i,β

Each element of the matrix can be written as: w i e i,γ e i,α e i,β

The last term is zero, so:

T c 2 c 2 s I = N i=0
f eq i c 2 t e i .e i

A.2 Advection-Diffusion Equation

A.2.1 1D Formulation

In 1D, the advection-diffusion process can be written as,

∂T ∂t + u ∂T ∂x = d ∂ 2 T
∂x 2 where u is the velocity and d the diffusion coefficient. The Lattice Boltzmann equation for advection-diffusion problem is still the same: f i (x + e i ∆x, t + ∆t) = f i (x, t) + w[f eq i (x, t) -f i (x, t)] (A.1) with:

f eq i (x, t) = ω i T (x, t)[1 + e i u cc 2 s ]
where c = ∆x ∆t

Chapman-Enskog Expansion

The advection-diffusion equation can be scaled spatially : x is set to x/ε, where ε is a small parameter:

∂T ∂t + uε ∂T ∂x = dε 2 ∂ 2 T ∂x 2 Now, if we scale time t with t/ε 2 as it is done for diffusion equation, then the first term on the left-hand side of the equation balances the diffusion term on the right hand side of the equation, hence the effect of advection term will be neglected. If we scale time t with 1/ε, in this case there will be a balance between the first and the second term of the left-hand side of the equation, and the effect of diffusion term will be neglected. To resolve this issue, we need to introduce dual time scales, i.e:

∂ ∂t = ε ∂ ∂t 1 + ε 2 ∂ ∂t 2
Hence, the first term in the above equation balances the advection term and the second term balances the diffusion term.

For one-dimensional problem with N lattice velocity components (c i , i = 1, 2, .., N ), the temperature can be expressed as,

T (x, t) = N i=1 f i (x, t)
Equation A.1 can be expanded using Taylor series:

f i (x+c i ∆t, t+∆t) = f i (x, t)+ ∂f i ∂t ∆t+ ∂f i ∂x c i ∆t+ ∆t 2 2 ( ∂ 2 f i ∂t 2 +2 ∂ 2 f i ∂t∂x c i + ∂ 2 f i ∂x 2 c 2 i )+0(∆t) 3
Introducing scaling, keeping term up to ε 2 we get:

f i (x + ce i ∆t, t + ∆t) = f i (x, t) + ε ∂f i ∂t 1 ∆t + ε 2 ∂f i ∂t 2 ∆t + ε ∂f i ∂x ce i ∆t + ∆t 2 2 (ε 2 ∂ 2 f i ∂t 2 1 + 2ε 2 ∂ 2 f i ∂t 1 ∂x ce i + ε 2 ∂ 2 f i ∂x 2 c 2 e 2 i ) + 0(∆t) 3
Substituting this equation in equation A.1:

ε ∂f i ∂t 1 ∆t+ε 2 ∂f i ∂t 2 ∆t+ε ∂f i ∂x ce i ∆t+ ∆t 2 2 (ε 2 ∂ 2 f i ∂t 2 1 +2ε 2 ∂ 2 f i ∂t 1 ∂x ce i +ε 2 ∂ 2 f i ∂x 2 c 2 e 2 i )+0(∆t) 3 = w[f eq i (x, t)-f i (x, t)]
The distribution function can be expanded in terms of the small parameter ε as,

f i (x, t) = f 0 i (x, t) + εf 1 i (x, t) + ε 2 f 2 i (x, t) + O(ε 3 )
Collecting the term of order ε 0 gives:

f 0 i (x, t) = f eq i (x, t) Since: N i=1
f eq i (x, t) = T (x, t)

Then:

N i=1 f 0 i (x, t) = T (x, t)
As:

N i=1 f i (x, t) = T (x, t)
This first equation gives:

N i=1 f j i (x, t) = 0 f or j = 0
Collecting the term of order ε 1 gives: since the right hand side is zero, and according to the momentum equations:

∂f 0 i ∂t 1 ∆t + ∂f 0 i ∂x ce i ∆t = -wf 1 i = -
N i=1 f 0 i = T N i=1 f 0 i c i = T u
The above equation formulate as:

∂T ∂t 1 + ∂(uT ) ∂x = 0 (A.3)
which is the advection equation. Latter we will need the derivative of f 1 i with respect to t 1 and x:

- 1 τ ∂f 1 i ∂t 1 = ∂ 2 f 0 i ∂t 2 1 + ∂ 2 f 0 i ∂x∂t 1 ce i
and:

- 1 τ ∂f 1 i ∂x = ∂ 2 f 0 i ∂t 1 ∂x + ∂ 2 f 0 i ∂x 2 ce i
These two equation gives:

- 1 τ ( ∂f 1 i ∂t 1 + ce i ∂f 1 i ∂x ) = ∂ 2 f 0 i ∂t 2 1 + 2 ∂ 2 f 0 i ∂t 1 ∂x ce i + ∂ 2 f 0 i ∂x 2 ce 2 i (A.4)
Collecting the term of order ε 2 gives:

∂f 1 i ∂t 1 + ∂f 0 i ∂t 2 + ∂f 1 i ∂x ce i + ∆t 2 ( ∂ 2 f 0 i ∂t 2 1 + 2 ∂ 2 f 0 i ∂t 1 ∂x ce i + ∂ 2 f 0 i ∂x 2 c 2 e 2 i ) = - 1 τ f 2 i Using equation A.4 ∂f 0 i ∂t 2 + ( ∆t 2 -τ )( ∂ 2 f 0 i ∂t 2 1 + 2 ∂ 2 f 0 i ∂t 1 ∂x ce i + ∂ 2 f 0 i ∂x 2 c 2 e 2 i ) = - 1 τ f 2 i
Hence:

∂f 0 i ∂t 2 + ( ∆t 2 -τ )( ∂ ∂t 1 ( ∂f 0 i ∂t 1 + ∂f 0 i ∂x c i ) + ∂ 2 f 0 i ∂t 1 ∂x ce i + ∂ 2 f 0 i ∂x 2 c 2 e 2 i ) = - 1 τ f 2 i
Summing the above equation over all states:

N i=1 ∂f 0 i ∂t 2 +( ∆t 2 -τ )( ∂ ∂t 1 ( N i=1 ∂f 0 i ∂t 1 + N i=1 ∂f 0 i ∂x ce i )+ N i=1 ∂ 2 f 0 i ∂t 1 ∂x ce i + 2 i=1 ∂ 2 f 0 i ∂x 2 c 2 e 2 i ) = - 1 τ 2 i=1 f 2 i
The right hand side is zero as well as the second term of the left hand side according to equation A.2. So by adding these 2 equations A.3 + A.5, we get:

∂T ∂t + ∂(uT ) ∂x + ε( ∆t 2 -τ ) ∂ 2 uT ∂t 1 ∂x = c 2 c 2 s (τ - ∆t 2 ) ∂ 2 T ∂x 2
There is an error term: ε( ∆t 2 -τ ) ∂ 2 uT ∂t 1 ∂x . The relaxation parameter can be related to the diffusion coefficient as:

d = (τ - ∆t 2 )c 2 s c 2

A.2.2 3D Formulation

In 3D, the advection-diffusion process can be written as,

∂T ∂t + u.∇T = ∇.( d∇T
)
where u is the velocity and d the diffusion coefficient. The Lattice Boltzmann equation for advection-diffusion problem is still the same: .6) where p = (x, y, z) and ∆p = ∆x = ∆y = ∆z with:

f i (p
f eq i (p, t) = ω i T (p, t)[1 + e i .u cc 2 s ]
where c = ∆p ∆t

Chapman-Enskog Expansion

The advection-diffusion equation can be scaled spatially : p is set to p/ε, where ε is a small parameter:

∂T ∂t + εu.∇T = ε 2 ∇.(d∇T )
Now, as in 1D, we introduce dual time scales, i.e:

∂ ∂t = ε ∂ ∂t 1 + ε 2 ∂ ∂t 2
For thre-dimensional problem with N lattice velocity components (c i , i = 0..N ), the temperature can be expressed as, T (x, t) = N i=0 f i (x, t) Equation A.6 can be expanded using Taylor series: The distribution function can be expanded in terms of the small parameter ε as:

f i (p + ce i ∆t, t + ∆t) = f i (p,
f i (p, t) = f 0 i (p, t) + εf 1 i (p, t) + ε 2 f 2 i (p, t) + O(ε 3 )
Collecting the term of order ε 0 gives: since the right hand side is zero, and according to the momentum equations:

f 0 i (p, t) = f eq i (p,
N i=0 f 0 i = T N i=0 f 0 i c i = T u
Indeed it is assumed here that the fluid is incompressible:

∇.u = 0

The above equation formulate as:

∂T ∂t 1 + u.∇T = 0 which is the advection equation. Latter we will need the derivative of f 1 i with respect to t 1 and space:

- 1 τ ∂f 1 i ∂t 1 = ∂ 2 f 0 i ∂t 2 1 + ∂ ∂t 1 (ce i .∇f 0 i ) and: - 1 τ ce i .∇f 1 i = ∂ ∂t 1 (cei.∇f 0 i ) + c 2 ∇∇ : f i t e i .e i )
These two equation gives:

- 1 τ ( ∂f 1 i ∂t 1 + ce i .∇f 1 i ) = ∂ 2 f 0 i ∂t 2 1 + 2 ∂ ∂t 1 (cei.∇f 0 i ) + c 2 ∇∇ : f i t e i .e i )
Collecting the term of order ε 2 gives:

∂f 1 i ∂t 1 + ∂f 0 i ∂t 2 + ce i .
f i t e i .e i )) = - 1 τ N i=0 f 2 i
The right hand side is zero as well as the second term of the left hand side according to equation A.7.

(uT )]

There is an error term: ε(τ -∆t 2 ) ∂ ∂t 1 [∇.(uT )] . The relaxation parameter can be related to the diffusion coefficient as:

d = (τ - ∆t 2 )c 2 s c 2
Let us consider the error term:

E = ε(τ - ∆t 2 ) ∂ ∂t 1 [∇.(uT )] = ε d c 2 s c 2 ∂ ∂t 1 [∇.(uT )]
To interpret this term, let us consider u as a constant, by taking the gradient of the advection equation we found and multiplying by u , we get:

εu.∇( ∂T ∂t 1 ) + ε 2 u 2 ∇ 2 T = 0
Thus the error term can be written as:

E = -ε 2 d c 2 s c 2 u 2 ∇ 2 T
And then the second equation becomes:

∂T ∂t + ∇.(uT ) = d(1 -u 2 c 2 s c 2 )∇ 2 T Therefore, the error term amounts to an undesired correction to the diffusion coefficient, proportional to the square of the advection velocity [Chopard 2009]. Parameters d and u are given by the physical problem, they are independent of the chosen time and space discretization ∆p and ∆t. Thus the correction to the diffusion coefficient scale as ( ∆t ∆p ) 2 . If we take the limits ∆t → 0 and ∆p → 0, with ∆p 2 /∆t constant (natural limit for a diffusion process), then the correction to the diffusion coefficient is of order O(∆t). The scheme is thus first order accurate in time and second order accurate in space. But if we take ∆t → 0 and ∆p → 0, with ∆p/∆t constant, then the error is of order O(ε 0 ) and the scheme is zeroth order accurate.

The error term can be eliminated if we consider:

Ω(f i ) = 1 τ (f eq i -f i ) + t i c i .δj i

A.3 Implementation

We use D3Q7 configuration, and for each lattice site, the seven distribution functions are calculated according to :

f i (r + c i ∆t, t + ∆t) = f i (r, t) + w[f eq i (r, t) -f i (r, t)] where S(T ) = R(T -T 0 ) in our case but can be any fonction of T.

In this case the lattice Boltzmann equation is changed to: f i (p + e i ∆p, t + ∆t) = f i (p, t) + w[f eq i (p, t) -f i (p, t)] + ∆tω i S (A.8)

Chapman-Enskog Expansion

The advection-diffusion-reaction equation can be scaled spatially : p is set to p/ε, where ε is a small parameter and a dual time scale is introduced:

ε ∂T ∂t 1 + ε 2 ∂T ∂t 2 + εu.∇T = ε 2 ∇.(d∇T ) + S
The fonction S should also be expanded in terms of the small parameter ε as:

S = S 0 + εS 1 + ε 2 S 2 + O(ε 3 )
Introducing this expression in the above equation yields:

ε ∂T ∂t 1 + ε 2 ∂T ∂t 2 + εu.∇T = ε 2 ∇.(d∇T ) + S 0 + εS 1 + ε 2 S 2
Clearly, it can be noted that S 0 = 0.

For thre-dimensional problem with N lattice velocity components (c i , i = 0..N ), the temperature can be expressed as:

T (x, t) = N i=0 f i (x, t)
Equation A.8 can be expanded using Taylor series: The distribution function can be expanded in terms of the small parameter ε as:

f i (p + ce i ∆t, t + ∆t) = f i (p,
f i (p, t) = f 0 i (p, t) + εf 1 i (p, t) + ε 2 f 2 i (p, t) + O(ε 3 )
Collecting the term of order ε 0 gives: Collecting the term of order ε 1 gives:

f 0 i (p, t) = f
∂f 0 i ∂t 1 ∆t + ∆tce i .∇f 0 i = -wf 1 i + ∆tω i S 1 = - ∆t τ f 1 i + ∆tω i S 1
Hence:

∂f 0 i ∂t 1 + ce i .∇f 0 i = - 1 τ f 1 i + ω i S 1
By taking the summation over all the values of i, we obtain: Finally:

N i=0 ∂f 0 i ∂t 1 + ∇( N i=0 ce i f 0 i ) = S 1 (A.9)
Since the first term of right hand side is zero, and according to the momentum equations:

N i=0 f 0 i = T N i=0 f 0 i c i = T u
Indeed it is assumed here that the fluid is incompressible:

∇.u = 0

The above equation formulate as:

∂T ∂t 1 + u.∇T = S 1 which is the advection equation. Latter we will need the derivative of f 1 i with respect to t 1 and space:

- 1 τ ∂f 1 i ∂t 1 = ∂ 2 f 0 i ∂t 2 1 + ∂ ∂t 1 (ce i .∇f 0 i )
and:

-

1 τ ce i .∇f 1 i = ∂ ∂t 1 (cei.∇f 0 i ) + c 2 ∇∇ : f i t e i .e i )
These two equation gives:

- 1 τ ( ∂f 1 i ∂t 1 + ce i .∇f 1 i ) = ∂ 2 f 0 i ∂t 2 1 + 2 ∂ ∂t 1 (cei.∇f 0 i ) + c 2 ∇∇ : f i t e i .e i ) (A.10)
Collecting the term of order ε 2 gives:

∂f 1 i ∂t 1 + ∂f 0 i ∂t 2 +ce i .∇f 1 i + ∆t 2 ( ∂ 2 f 0 i ∂t 2 1 +2 ∂ ∂t 1 (ce i .∇f 0 i )+c 2 ∇∇ : f 0 i t e i .e i ) = - 1 τ f 2 i +ω i S 2
Using equation A.10:

∂f 0 i ∂t 2 + ( ∆t 2 -τ )( ∂ 2 f 0 i ∂t 2 1 + 2 ∂ ∂t 1 (ce i .∇f 0 i ) + c 2 ∇∇ : f 0 i t e i .e i ) = - 1 τ f 2 i + ω i S 2
Hence:

∂f 0 i ∂t 2 +( ∆t 2 -τ )( ∂ ∂t 1 ( ∂f 0 i ∂t 1 +ce i .∇f 0 i )+ ∂ ∂t 1 (ce i .∇f 0 i )+c 2 ∇∇ : f 0 i t e i .e i ) = - 1 τ f 2 i +ω i S 2
Summing the above equation over all states:

N i=0 ∂f 0 i ∂t 2 + ( ∆t 2 -τ )( ∂ ∂t 1 ( N i=0 ∂f 0 i ∂t 1 + c N i=0 e i .∇f 0 i ) + ∂ ∂t1 ( N i=0 ce i .∇f 0 i ) + c 2 ∇∇ : N i=0 f i t e i .e i )) = - 1 τ N i=0 f 2 i + S 2
The first term of the right hand side is zero and according to equation A.9: where e ᾱ = -e α and fi (r, t) = f i (r, t) + w[f eq i (r, t) -f i (r, t)]

N i=0
Computational Modeling of Radiofrequency Ablation for the Planning and Guidance of Abdominal Tumor Treatment Abstract: The outcome of radiofrequency ablation (RFA) treatment for abdominal tumors is challenged by lack of patient-specific planning. In particular, the presence of blood vessels, time-and spatial-varying thermal conductivity makes the prediction of the extent of ablated tissue difficult. This results in possible incomplete treatments and increased risks of recurrence and may compromise the overall outcome. By providing predictive tools, biophysical models can help clinicians to plan and guide the procedure for an effective treatment. For such models to be useful in clinical routine, the computation has to be fast and accurate. We first present a detailed computational model of the biophysical mechanisms involved in RFA of hepatic tumors such as heat diffusion, cellular necrosis and hepatic blood flow. This model simulates the extent of ablated tissue given the probe location and some biological parameters such as tissue conductivity or heat capacity. The model is based on medical images, from which patient-specific models of the liver, visible vessels and tumors are build. A new approach for solving these partial differential equations based the Lattice Boltzmann Method is introduced. The model is implemented with a material speed-up thanks to the graphics card to reach near real-time computation. The resulting simulations are thoroughly verified against an analytical solution.

In a second step, we aim at validating the model. It is first evaluated against clinical data of patients who underwent RFA of liver tumors. The importance of taking into account the blood flow and of estimating the parameters is highlighted. Then, a comprehensive pre-clinical experiment combining multi-modal, pre-and post-operative anatomical and functional images, as well as the interventional monitoring of the temperature and delivered power is presented. This enables an end-toend pre-clinical validation framework that considers the most comprehensive data set for model validation.

Then, we automatically estimate patient-specific parameters using inverse problem algorithms, to better predict the ablated tissue. This personalization strategy has been tested and evaluated on seven ablations from three clinical cases. From the pre-clinical study, we can go further in the personalization by comparing the simulated temperature and delivered power with the actual temperature and delivered power measured during the procedure. This second personalization strategy has been tested and evaluated on twelve ablations from five pre-clinical cases. These contributions have led to promising results, and open new perspectives in RFA guidance and planning.

4. 7

 7 Two hypotheses for the probe position on patient 06 . . . . . . . . . 60 4.8 Evaluation of the effect of the RFA probe. . . . . . . . . . . . . . . . 60 4.9 Effect of the probe position on patient 06 . . . . . . . . . . . . . . . 61 4.10 Evaluation of the effect of the segmentation. . . . . . . . . . . . . . . 62 4.11 Two segmentation results of patient 03 . . . . . . . . . . . . . . . . . 62 4.12 Effect of the segmentation on RFA simulation for patient 03 . . . . . 62 4.13 Effect of the diffusivity parameter on the error for patient 02 . . . . 63 4.14 Parameter tuning on patient 02 . . . . . . . . . . . . . . . . . . . . . 64 5.1 Parameter estimation pipeline. . . . . . . . . . . . . . . . . . . . . . 71 5.2 Synthetic case for CFD solver validation. . . . . . . . . . . . . . . . . 72 5.3 Synthetic set-up for the verification of the optimization framework. . 73 5.4 Convergence analysis of the optimization framework. . . . . . . . . . 75 5.5 Results of the parameter estimation on the computed lesion. . . . . . 76 5.6 Box plot of the error on the 7 ablations with personalized parameters. 77 5.7 Parameters personalized on three patients. . . . . . . . . . . . . . . . 77 5.8 Comparison of personalized computation with post-operative registered ground truth for patient 1-1 . . . . . . . . . . . . . . . . . . . . 78 5.9 Comparison of personalized computation with post-operative registered ground truth for patient 1-2 . . . . . . . . . . . . . . . . . . . . 79 5.10 Comparison of personalized computation with post-operative registered ground truth for patient 1-3 . . . . . . . . . . . . . . . . . . . . 80 6.1 Pre-operative CT and Post-operative MRI available for Pig 5. . . . . 82 6.2 Pipeline presenting the available pre-clinical data. . . . . . . . . . . . 83 6.3 Surrogate tumor implantation. . . . . . . . . . . . . . . . . . . . . . 84 6.4 RFA probe inside the liver and on a U.S image. . . . . . . . . . . . . 85 6.5 CFD Model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 6.6 Cell state evolution with the adapted cellular necrosis model. . . . . 88 6.7 Comparison of the 12 delivered power curves. . . . . . . . . . . . . . 89 6.8 Delivered power over ablation : the RFA protocol is not exactly the same for all pigs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 6.9 Probe: mesh model used and interventional MRI data. . . . . . . . . 91 6.10 Segmentations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 6.11 Smooth vessel segmentation. . . . . . . . . . . . . . . . . . . . . . . . 93 6.12 Velocity from Phase-Contrast MRI. . . . . . . . . . . . . . . . . . . . 94 6.13 Porosity map used for Pig 1. . . . . . . . . . . . . . . . . . . . . . . . 95 6.14 Evaluation of the Post-to-Pre registration. . . . . . . . . . . . . . . . 96 7.1 Results of the Post-to-Pre registration on the liver of Pig 1 liver. . . 101 7.2 Post-to-Pre registration of arterial mesh of Pig 2. . . . . . . . . . . . 102 7.3 Post-to-Pre registration of liver mesh of Pig 2. . . . . . . . . . . . . . 103 7.4 Evaluation of the vessels smoothing of Pig 4 . . . . . . . . . . . . . . 104 7.5 Synthetic case for the parameter estimation. . . . . . . . . . . . . . . 105 7.6 Computed delivered power and temperature on the synthetic case. . 7.7 Convergence curves of 2 optimization algorithms. . . . . . . . . . . . 7.8 Sensitivity analysis on the heat capacity of undamaged cells c U t . . . . 7.9 Sensitivity analysis on the heat capacity of vulnerable cells c V t . . . . 7.10 Sensitivity analysis on the heat capacity of necrotic cells c N t . . . . . . 7.11 Sensitivity analysis on the heat conductivity dt . . . . . . . . . . . . . 7.12 Evaluation on synthetic data. . . . . . . . . . . . . . . . . . . . . . . 7.13 Computed delivered power and temperature for tumor 4-3. . . . . . . 7.14 Zoom on the cooling phase of the temperature for tumor 4-3. . . . . 7.15 Simulated necrosis versus ground truth for tumor 4-2. . . . . . . . . 7.16 Box plot of the error on the 12 ablations from 5 pigs. . . . . . . . . . 7.17 Box plot of the similarity scores computed on the volume of the 12 ablations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.18 Box plot of the error in term of delivered power and temperature on the 12 ablations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.19 Simulated delivered power and temperature for tumor 3-1 with nominal parameters and with personalized parameters. . . . . . . . . . . 7.20 Simulated delivered power and temperature for tumor 3-1 with and without the advection term. . . . . . . . . . . . . . . . . . . . . . . . 7.21 Parameters personalized on two tumors from two different pigs. . . . 7.22 Personalization of the delivered power and temperature curves for tumor 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.23 Zoom on the cooling temperature curves for tumor 1. . . . . . . . . . 7.24 Effect of the advection on the simulated necrosis versus ground truth for Pig 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.25 Effect of the location of the ablation on the temperature dispersion. . 7.26 Effect of the advection on the delivered power and the temperature curves for tumor 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.27 Zoom on the cooling temperature curves for tumor 1 computed with the advection. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.28 Effect of the registration on the simulated lesion on Pig 2. . . . . . . 7.29 Simulated lesion versus ground truth for tumor 3-2. . . . . . . . . . . 7.30 Cell death model versus thermal dose computation. . . . . . . . . . . 7.31 The simulated necrotic extent around tumor 3-3 is subcapsular. . . . 7.32 MR thermometry acquired during the ablation of tumor 4-1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.

Figure 1

 1 Figure 1.1: A schematic diagram showing how the chapters are connected and the different available data used in each chapter.

Figure 2 . 1 :

 21 Figure 2.1: Anatomy of the human liver. Image from the internet 1 .

Figure 2

 2 Figure 2.2: A schematic diagram showing an RF needle deployed in a liver tumor.A complete electrical circuit is then created through the patient body allowing the alternating current of RF energy to go between the needle active electrode and the grounding pads. The white area is the targeted tumor, and the surrounding orange area, the necrosis lesion produced by RFA. Image from[START_REF] Curley 2003steven | Radiofrequency ablation of malignant liver tumors[END_REF]].

Figure 2

 2 Figure 2.3: Drawing depicts how the heat expands around the needle electrodes into the surrounding tissue to form a thermal sphere during RFA intervention. In this case, the RFA probe just has 4 tips. Image from Internet 4 .

Figure 2 . 4 :

 24 Figure 2.4: Liver extraction showing RFA-induced lesion after in vivo ablation on a pig liver.

Figure 2 . 5 :

 25 Figure 2.5: A CT slice showing radiofrequency ablation probe in the liver. Image from Wikipedia 5

Figure 2 . 7 :

 27 Figure 2.7: RITA probe representation from [Payne 2011].

Figure 3

 3 Figure 3.1: Steps of the proposed method (green: input, blue: processes, red: output).

  Figure 3.3: Model of the hepatic circulatory system. Arrows denote blood flow. Circles and squares denote portal and hepatic vessel tips. See text for details.

FigureFigure 3 . 5 :

 35 Figure 3.5: Set-up and results of the hepatic blood flow computation with a zoom inside the vena cava on the right and at the extremities of the hepatic veins on the left.

  .6) and Neumann boundary conditions is employed.

Figure 3 . 6 :

 36 Figure 3.6: D3Q19 scheme used in the LBM computation of the blood flow.

  Fig 3.7 shows the results of the two different computations on a patient-specific geometry.

Figure 3 . 7 :

 37 Figure 3.7: (Left): Blood flow computation with the Decoupled FEM Porous-Vessels CFD Model; (Right): Blood flow computation with the Fully Coupled LBM CFD Model.In both cases, only the flow going through the portal inlet circled in white are shown. We can see differences between the results given by both computations. The area circled in black shows flow lines, which appear with the first computation and not the second one.

Figure 3

 3 Figure 3.8: Cell state evolution over time when tissue is heated at 105 • C during 10 minutes.

Figure 3 .

 3 Figure 3.10: Layout of the isotropic Cartesian grid for the implementation of LBM boundary condition with a curved wall boundary.

Figure 3

 3 Figure 3.11: Scheme showing the modeling of the RFA probe with a sphere centered at the center of the tumor.

  Figure 3.12: A detailed anatomical model of the liver is estimated from a standard clinical CT image. See text for details.

Figure 3

 3 Figure 3.13: Synthetic set-up for the verification against analytical solution: (Left): Synthetic set-up used in the quantitative verification of the framework at time t = 0 s. (Right): Initial temperature distribution along the x-axis of the cuboid domain. See text for details.

Figure 3 .

 3 Figure 3.14: (Left): Spatial convergence analysis for a fixed time step of 0.01 s for one point in the domain described in Figure 3.13. As one can see, the proposed framework quickly converges to the analytical solution of the advection-diffusion equation. (Right): The computed solution for a resolution of 1 mm and a time-step of 0.075 s compared to the analytical solution for the same point. The RMS error is reduced to 0.26 • C.

  3.1).

  Figure 3.15: (Left): Computation time for simple CPU implementation without any parallelization, for multithreading implementation with OpenMP and 22 threads and GPU implementation with CUDA. (Right): Computation speed-up with respect to the number of threads used (parallel optimization with OpenMP).

Figure 4 .

 4 Figure 4.1: (Left) Pre-operative image of patient 04. The tumor is visible. (Right) Post-operative image for patient 04, the necrotic area is visible around hepatic vein.

Figure 4 .

 4 Figure 4.2:The registration is performed using masks with vessels and parenchyma as landmarks. The post-operative mask is registered (Right) to the pre-operative mask (Left), with an elastic registration algorithm using ANTS[Avants 2009]. Presented here are the results for patient 10.

4. 4 Figure 4 . 3 :

 443 Figure 4.3: Results of the computation for patient 04, the streamlines represent the parenchyma flow and are color-coded with respect to the velocity magnitude. As one can see the lesion computed with advection follows the blood flowing from the portal vein to the sushepatic veins. See text for details.

Figure 4 .

 4 3 and Figure 4.5). Cell death area computed using the model compared qualitatively well with the observed post-operative necrosis zone as we can see on Figure 4.3 and Figure 4.5. Quantitatively, in most of the cases, average point-to-mesh errors were within clinical acceptance as they were significantly lower than the different diameter configurations of the probes (3 to 5 cm), as shown on Figure 4.4. The values are reported inTable 4.1.

Figure 4 . 4 :Figure 4 . 5 :

 4445 Figure 4.4: Box plot of the error on the 14 ablations. (Left): Point-to-Mesh error in mm computed on the surface of the lesion. (Right): Similarity scores in % computed on the volume of the lesion.

Figure 4 . 6 :

 46 Figure 4.6: Box plot of the error on the 14 ablations computed without the advection term. (Left): Point-to-Mesh error in mm computed on the surface of the lesion. (Right): Similarity scores in % computed on the volume of the lesion.

Figure 4 . 7 ,Figure 4 . 7 :

 4747 Figure 4.7: The 2 different hypotheses used for the probe placement (patient 06). (Hypothesis 1 ): The probe is placed at the center of the tumor = regular hypothesis. (Hypothesis 2 ): The probe is placed at the center of the post-op registered lesion.

Figure 4 . 8 :Figure 4 . 9 :

 4849 Figure 4.8: Evaluation of the effect of the probe placement with the two hypotheses described in Figure 4.7.

Figure 4 .Figure 4 .

 44 Figure 4.10: Evaluation of the effect of the segmentation on the simulated lesion for patient 03.

Figure 4 .

 4 Figure 4.13: Effect of the diffusivity parameter on the mean and standard deviation of the point to mesh error, the DICE, Sensitivity and PPV scores (patient 02). D x 8 gives the larger DICE score and the smaller point-to-mesh error.

Figure 4 .

 4 Figure 4.14: (Left): computed lesion in red is smaller than the post-op lesion in blue. (Right): A better match is achieved between the computed lesion and the ground truth with adjusted heat conductivity for the same patient (patient 02).

Figure 5

 5 Figure 5.1: Estimation of the personalized parameters and the forward model (blue: input, green: processes, purple: output).

Figure 5 . 2 :

 52 Figure 5.2: Synthetic case used to validate the CFD solver. (Top Left): A cylindrical mesh with spherical inlet and outlet is used. The boundary conditions are the inflow and the outlet pressure. (Top Right): The velocity field given by the CFD solver in the second case. (Down Left): The porosity field used in the first case. (Down Right): The porosity field used in the second case.

Figure 5

 5 Figure 5.3: Set-up of the synthetic case for the verification of the optimization framework. (Top): The cylinder with the porosity field used. The boundary conditions are the output pressure and the input flow. (Down): The heat distribution initially applied and the velocity distribution used.

Figure 5 . 4 :

 54 Figure 5.4: Convergence analysis of the optimisation framework. We use DAKOTA and manage to estimate two parameters ( dt and ε) in 32 minutes after 36 iterations with a mean of the point-to-mesh error of 10 -3 mm.

Figure 5 .Figure 5 . 5 :

 555 Figure 5.5: Computed necrosis compared qualitatively well with the predicted lesion after personalization on the first tumor of each patient.

Figure 5 . 6 :Figure 5 . 7 :

 5657 Figure 5.6: Box plot of the error on the 7 ablations with personalized parameters. (Left): Point-to-Mesh error in mm computed on the surface of the lesion). (Right): DICE score in % computed on the volume of the lesion.

Figure 5 4 : 5 :Figure 5 . 9 :

 54559 Figure 5.8: Parameter Estimation on the first tumor of the first patient (patient 1-1). (Left): The post-operative lesion appears around the tumor. (Middle-Left): Computed necrosis with personnalised parameters. (Middle-Right): Overlay of the computed lesion with personnalized parameters and ground truth. (Right): Overlay of the computed lesions with nominal parameters, with personalized parameters and ground truth. In each case, the blood flow in the area of interest is represented by streamlines in white.

Figure 6 .

 6 Figure 6.1: (Left) Venous phase of the pre-operative CT of Pig 5. The surrogate tumor is visible, as well as the gall bladder. (Right) Post-operative MRI T2 for Pig 5, the necrotic area is visible around the surrogate tumor, as well as the gall bladder.

Figure 6 . 2 :

 62 Figure 6.2: Pipeline presenting the available clinical data in blue, the pre-processing needed to use those data (segmentation in green, registration in red), the different computations done in grey and finally a possible personalization framework.

Figure 6

 6 Figure 6.3: Two CT images of Pig 1 before and after surrogate tumors implantation. (Left): Before tumor implantation; (Right): After tumor implantation.

Figure 6 . 4 :

 64 Figure 6.4: (Left): Photo of the probe inside the pig liver; (Right): U.S image of the probe next to the surrogate tumor. The probe diameter is shown in red, the tumor diameter in green.

Figure 6 . 5 :

 65 Figure 6.5: Model of the hepatic circulatory system. Black arrows (resp. grey arrows) denote blood flow (resp. outlet pressures) which are fixed by the boundary conditions. See text for details.

Figure 6 . 6 :

 66 Figure 6.6: Cell state evolution over time when tissue is heated at 105 • C during 10 minutes with the adapted cellular necrosis model proposed in this chapter.

Figure 6 . 7 :

 67 Figure 6.7: Plot of the power delivered for the 12 ablations for the first 200 seconds. The mean and standard deviation are computed at time t = 50 s.

Figure 6 .

 6 Figure 6.7 represents the power delivered for the 12 ablations for the first 200 seconds. However, the applied RFA protocol is not exactly the same for all ablations. Eight ablations are performed through several short cooling and heating periods, in this case, the delivered power is reduced during the different cooling stages. On the other hand, the other four ablations included only one long final cooling stage after a continuous heating period (red, blue, cyan and green curves in Figure. 6.7). Figure.6.8 shows the four cases whose the delivered power over ablation is a long heating phase followed by a cooling stage, and one typical case whose ablation is performed through several short cooling and heating periods (pig 5 -tumor 1).

Figure 6 .

 6 Figure 6.8: (Left): Delivered power for the four ablations included only one long final cooling stage after a continuous heating period; (Right): Delivered power for pig5 -tumor 1, for which ablation was performed through several short cooling and heating periods. In both case, the grey area represents the domain considered in Figure. 6.7

Figure 6 . 9 :

 69 Figure 6.9: (Left): Mesh model of the probe with the 9 tips derived from a CT image of the probe only; (Right): Anatomical MRI with the MRI-compatible probe implanted in the liver. Four tips of the probe are visible.

Figure 6 .

 6 Figure 6.10: Segmentation of the liver in magenta, arterial vessels in red, portal and hepatic networks in light blue and dark blue respectively, surrogate tumors in dark green and gall bladder in light green. (Left): Overlay of the segmented areas on pre-operative CT image. (Middle): Creation of the 3D model. (Right): Vessels and liver parenchyma.

Figure 6 .

 6 Figure 6.11: (Left): Centerlines extraction from the vessel segmentation. The segmentation gives non-smooth vessels due to the anisotropy and the low resolution of the input image. (Right): Smooth vessels generated from the centerlines.

Figure 6 .

 6 Figure 6.14: (Left): Superposition of the post-operative mesh model in red on the pre-operative binary image; (Right): Comparison between the registered postoperative vessel mesh and the pre-operative one.

Figure 7 . 5 :Figure 7 . 6 :

 7576 Figure 7.5: Set-up of the synthetic case. (Left): The temperature map around one tip of the probe. (Right): In red, the necrotic area created around the probe in dark.

Figure 7 . 7 :

 77 Figure 7.7: Convergence curves of the COBYLA algorithm in blue and of the pattern search algorithm in green. The two optimization algorithms are compared on a synthetic case.

Figure 7

 7 Figure 7.8: Sensitivity analysis on the heat capacity of undamaged cells c U t . It affects strongly the time of the peak value of the delivered power. (Left): Effect on the delivered power. (Right): Effect on the cooling temperature.

Figure 7 . 9 :

 79 Figure 7.9: Sensitivity analysis on the heat capacity of vulnerable cells c V t . (Left): Effect on the delivered power. (Right): Effect on the cooling temperature.

Figure

  Figure 7.10 illustrates the effect of the heat capacity of necrotic cells, c N t on the delivered power and on the cooling temperature distribution. This parameter does not largely affect the delivered power neither the cooling temperature.Finally, Figure7.11 illustrates the effect of the heat conductivity, dt on the delivered power and on the cooling temperature distribution. The heat conductivity has a clear impact on the cooling temperature distribution.

Figure 7 .

 7 Figure 7.10: Sensitivity analysis on the heat capacity of necrotic cells c N t . (Left): Effect on the delivered power. (Right): Effect on the cooling temperature.

Figure 7 .

 7 Figure 7.11: Sensitivity analysis on the heat conductivity dt . (Left): Effect on the delivered power. (Right): Effect on the cooling temperature.

Figure 7 .

 7 Figure 7.13: (Left): The computed delivered power curves over ablation are compared with the measured one. (Right): The temperature distribution after ablation of tumor 4-3. During the heating phase, five temperatures are directly imposed from the five thermistors measurements (superposition of green and black curves)and the four remaining temperatures are imposed from a linear interpolation from those measurements (four green curves). After the use of personalized parameters, the error between the measurement and the computation is reduced from 18.9 to 8.8 W for the power and from 3.5 to 3.0 • C for the temperature. A zoom on the cooling phase is presented in the following Figure7.14.

Figure 7

 7

Figure 7 .

 7 Figure 7.15: Comparison between the simulated and the post-operative necrotic areas on tumor 4-2. (Left): The simulated lesion is showed around the RFA probe in the subject-specific geometry. (Right): Zoom on the ablation area, the simulated lesion is qualitatively close to the registered post-operative lesion. For this tumor, the point-to-mesh error is 4.6 ± 3.9 mm.

Figure 7 .

 7 Figure 7.16: Point-to-Mesh error in mm computed on the surface of the lesion for three different computations: with nominal parameters, with personalized parameters and with the advection term. For the computation with the advection term, only 6 ablations are considered. (Left): Mean of the point to mesh (Right): Standard deviation of the point to mesh error.

Figure 7 .Figure 7 .Figure 7 .

 777 Figure 7.17: Similarity scores in % computed on the volume of the lesion for three different computations: with nominal parameters, with personalized parameters and with the advection term. For the computation with the advection term, only 6 ablations are considered. (Left): DICE scores (Middle): Sensitivity (Right): Positive Predictive Value.

Figure 7 .

 7 Figure 7.20: Comparison between the simulated and the measured temperature on tumor 3-1 for 2 different computations: with and without the advection term. During the heating phase, the temperature is imposed, the power is simulated and compared to the measured one. The error in term of delivered power is constant (35.0 W versus 35.1 W). During the cooling phase, the simulated (non-imposed) temperature is compared to the measured one, whenever measurement is available (the RFA power was turned off during the cooling phases). There is no real difference in term of cooling temperature error: 6.3 • C versus 5.9 • C.

Figure 7 .

 7 Figure 7.21: Comparison of the parameters (the vulnerable rate coefficient of the cell death model δ, the heat capacity of undamaged cells c U t , the heat conductivity dt ) personalized on two tumors from two different pigs against the nominal values from the literature.

Figure 7 .

 7 Figure 7.23: Zoom on the cooling temperature distribution after ablation from tumor 1. The personalized parameters do not affect the cooling temperatures.

Figure 7 .

 7 Figure 7.24: Comparison between the simulated and the post-operative necrotic areas on Pig 1. (Left): The post-operative lesion is showed around the RFA probe in the subject-specific geometry. (Right): Zoom on the ablation area, the simulated lesions with and without advection are compared to the registered post-operative lesion.The advection has an effect on the shape of the lesion, however in this case, it does not improve the point-to-mesh error: 6.9 ± 5.2 mm versus 4.8 ± 4.4 mm.

Figure 7 .

 7 Figure 7.25: Effect of the location of the ablation on the temperature dispersion.(Left): Pig 1, the probe is close to vessels, the tips do not cool down similarly (Figure7.27). (Right): Tumor 3-1, the probe is far from vessels, the tips cool down similarly (Figure7.20). In this case, we could not segment any post-op lesion.

Figure 7 .

 7 Figure 7.26: Effect of the advection on (Left) the delivered power curve over ablation (Right) the temperature distribution after ablation from tumor 1. The error between the measurement and the computation is reduced from 10.3 • C to 9.7 • Cfor the cooling temperature with the advection, whereas it does not really affect the delivered power (19.9 W to 20.8 W). A zoom on the cooling phase is presented on the following Figure7.27.

Figure 7 .

 7 Figure 7.28: Effect of the registration on Pig 2 (Left): The post-operative lesion after the non-rigid registration is showed in the subject-specific geometry. The point-tomesh error is 8.4 ± 5.0 mm (Right): The post-op lesion after a single translation.The poinr-to-mesh error reduces to 3.6 ± 2.9 mm.

Figure 7 .

 7 Figure 7.29: Comparison between the simulated and the post-operative necrotic areas on tumor 3-2. (Left): The post-operative lesion is showed around 2 possible positions of the RFA probe. (Right): Zoom on the ablation area, the simulated lesion is closer to the registered post-operative lesion when the probe is in position 2. The point-to-mesh error is 12.6 ± 7.7 mm in Position 1 and 6.2 ± 4.0 mm in Position 2.

Figure 7 .

 7 Figure 7.30: Effect of the cell death model (Left): The simulated lesion computed using the cell death model is showed in the subject-specific geometry of Pig 2. (Right): The area where at least a thermal dose of 240 in equivalent minutes at 43 • C is shown. The post-op lesion appears in dark in both cases.

Figure 7 .

 7 Figure 7.31: (Left): The simulated necrotic extent around the RFA probe, close to the border of the liver from tumor 3-3. The post-operative lesion is shown in brown. (Right): The liver extraction confirms that the necrotic lesion is subcapsular.

Figure 7 .

 7 Figure 7.32: MR thermometric images acquired during the ablation of tumor 4-1. (Left): image acquired before the ablation, showing an isotropic and constant temperature around the RFA probe (liver temperature). (Right): image acquired during the ablation, when the RF power is decreased. It shows the temperature distribution around the RFA probe.

ω

  u] And let us consider:c i = ce iwhere c = dr/dt and e i are vectors with coefficient 1 or -1.We have the following properties (called the isotropic properties) which should be true to keep the lattices symmetries[Chai 2013]: i e T i .e i = c 2 s I Due to the last isotropic property, we have:

  i .e i = T N i=0 w i [1 + c i .u] t e i .e iAs:

  .e i ) α,β = T N i=0 w i [1 + c i .u]e i,α e i,β = T c 2 s δ α,β + T c

i

  By taking the summation over all the values of i, we obtain:

i

  (p, t) = 0 f or j = 0Collecting the term of order ε 1 gives: By taking the summation over all the values of i, we obtain:

  So by adding the first equation times ε and the second one times ε 2 , we get:∂T ∂t + ∇.(uT ) = c 2 c 2 s (τ -∆t 2 )∇ 2 T + ε(τ -

  -Diffusion-Reaction EquationWe are interested in the 3D advection-diffusion-reaction equation:∂T ∂t + u.∇T = ∇.(d∇T ) + S(T )

  i (p, t) = T (p, t) t) = 0 f or j = 0

  1.1 Clinical Context . . . . . . . . . . . . . . . . . . . . . . . . .

1 1.1.2 Role of Medical Imaging . . . . . . . . . . . . . . . . . . . . . 2 1.1.3 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 Main Contributions and Manuscript Organization . . . . . 3 1.2.1 Main Contributions . . . . . . . . . . . . . . . . . . . . . . . 3 1.2.2 Organization of the Thesis . . . . . . . . . . . . . . . . . . . . 4

  State of the Art of Existing Radiofrequency Ablation Models 22 2.5.1 Computational Models of RFA . . . . . . . . . . . . . . . . . 22 2.5.2 Model Simplification . . . . . . . . . . . . . . . . . . . . . . . 23 2.5.3 Model Validation . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.5.4 Model Personalization . . . . . . . . . . . . . . . . . . . . . .

	Chapter 2
	State of the Art on
	Radiofrequency Ablation of
	Hepatic Tumors
		25
	2.5.5 Surrogate Model for Optimization of the Probe Placement . .	26
	2.6 MR Thermometry . . . . . . . . . . . . . . . . . . . . . . . . . 27
	2.6.1 Thermal Tissue Properties Acquisition during HIFU based on
	MR Thermometry . . . . . . . . . . . . . . . . . . . . . . . .	27
	2.6.2 MR Thermometry Filtered by Bioheat Equation . . . . . . .	27

Contents 2.1 Structure and Function of the Liver . . . . . . . . . . . . . . 10 2.2 Radiofrequency Ablation . . . . . . . . . . . . . . . . . . . . . 11 2.2.1 Radiofrequency Ablation Procedure . . . . . . . . . . . . . . 11 2.2.2 Monitoring of RFA based on MR Thermometry . . . . . . . . 13 2.3 Percutaneous Tumor Ablation Techniques . . . . . . . . . . 14 2.3.1 High Intensity Focused Ultrasound . . . . . . . . . . . . . . . 15 2.3.2 Microwave Ablation . . . . . . . . . . . . . . . . . . . . . . . 15 2.3.3 Cryotherapy . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.3.4 Laser Ablation . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3.5 Electroporation . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.4 Technical Background . . . . . . . . . . . . . . . . . . . . . . . 17 2.4.1 Models of Electrical Heating . . . . . . . . . . . . . . . . . . 17 2.4.2 Biophysical Models of Heat Transfer in Tissue . . . . . . . . 18 2.4.3 Models of Cellular Necrosis . . . . . . . . . . . . . . . . . . . 20 2.5

Table 2 .

 2 1: Values from literature[Payne 2011] of the parameters used in the computation of RFA.

	parameter description	value
	ρ	blood and tissue densities	1.06 × 10 3 kg m -3
	c bl	blood heat capacity	4.18 × 10 3 J(kg K) -1
	c U ti , c V ti	tissue heat capacity	3.6 × 10 3 J(kg K) -1
	c N ti	tissue heat capacity in necrotic cells	0.67 × 10 3 J(kg K) -1
	d	blood and tissue heat conductivities	0.512 × (1 + 0.00161 × (T -310)) W(m K) -1
	H	convective transfer coefficient	24.4 × 10 5 W (m 3 K) -1
	ε	blood volume fraction	0.1
	κ	permeability	4.0 × 10 -11 m 2
	µ	dynamic viscosity of the blood	0.0035 P a s
	ϕ vcin	vena cava inflow	2.0 L min -1
	ϕ i	flow through the inlets of the hepatic veins	1.6 L min -1
	p 0	vena cava outlet pressure	3 mmHg
	kf	forward rate constant	3.33 × 10 -3 s -1
	k b	backward rate constant	7.77 × 10 -3 s -1
	T k	parameter of cell state model	40.5 • C

Table 4 .

 4 1: Evaluation of the RFA modeling on patient data (PPV: Positive Predictive Value).

	Patient	Probe diameter	Point-to-mesh	DICE	Sensitivity	PPV
		(cm)	error (mm)	(%)	(%)	(%)
	01	3	8.24 ± 8.48	54.3	46.0	66.2
	02	4 then 5	8.25 ± 5.92	61.7	45.6	95.6
	03	4	21.65 ± 16.96	17.7	22.7	14.5
	04	4 then 5	10.91 ± 8.67	31.4	89.5	19.0
	05-1	4	7.28 ± 7.03	45.0	80.8	31.2
	05-2	3.5	6.80 ± 6.25	54.2	44.8	68.5
	06	4 then 5	10.77 ± 5.01	37.6	90.7	23.7
	07-1	4 then 5	8.34 ± 6.79	40.8	66.4	29.4
	07-2	4 then 5	12.92 ± 12.32	45.2	63.8	35.0
	07-3	4 then 5	11.61 ± 11.92	40.3	93.6	25.7
	08-1	3	7.97 ± 5.41	30.6	32.1	29.1
	08-2	3	5.32 ± 4.76	61.6	65.5	58.2
	09	3	11.96 ± 8.50	21.3	98.9	11.9
	10	4 then 5	10.33 ± 11.20	44.1	96.8	28.5
	Mean	-	10.17 ± 8.52	41.8	66.9	38.3

Table 4 .

 4 2: Evaluation of the RFA modeling on patient data without advection.

	Patient	Point-to-mesh	DICE (%)	Sensitivity (%)	PPV (%)
		error (mm)
	01	9.03 ± 9.74	52.4	43.7	65.3
	02	8.00 ± 5.01	64.6	48.9	95.2
	03	24.44 ± 20.91	18.3	22.4	15.5
	04	10.57 ± 5.59	28.8	96.6	16.9
	05-1	5.56 ± 5.20	51.7	83.8	37.4
	05-2	6.19 ± 4.56	54.5	59.8	50.1
	06	8.25 ± 5.92	37.9	82.6	24.6
	07-1	6.06 ± 3.99	55.3	79.7	42.3
	07-2	7.98 ± 5.05	56.8	71.2	47.2
	07-3	4.73 ± 3.30	67.8	87.9	55.2
	08-1	6.70 ± 4.43	38.8	51.6	31.0
	08-2	5.07 ± 4.45	60.5	61.3	59.7
	09	6.91 ± 4.65	33.0	97.9	19.9
	10	6.76 ± 4.08	51.5	98.8	34.8
	Mean	8.30 ± 6.21	48.0	70.4	42.5

Table 4 .

 4

	Patient 06	Point-to-mesh error	DICE	Sensitivity	PPV
	Hypothesis 1	10.77 ± 5.01 mm	37.6 %	90.7 %	23.7 %
	Hypothesis 2	8.01 ± 6.50 mm	50.0 %	96.3 %	33.8%
		Point-to-mesh error
		(mean)
		(mm)			(%)

3: Evaluation of the effect of the probe placement with the two hypotheses described in Figure

4

.7.

Table 4 . 4 :

 44 Evaluation of the effect of the segmentation on the simulated lesion for patient 03.

	Patient 03	Point-to-mesh error	DICE	Sensitivity	PPV
	Segmentation 1	21.65 ± 16.96 mm	17.7 %	22.7 %	14.5 %
	Segmentation 2	19.07 ± 17.11 mm	22.8 %	28.4 %	19.0%

Table 4 . 5 :

 45 Estimation of the heat conductivity on patient 02 (PPV: Positive Predictive Value). D x 8 gives the larger DICE score and the smaller point-to-mesh error.

	Patient 02	Point-to-mesh	DICE (%)	Sensitivity (%)	PPV (%)
	Probe diameter:	error (mm)
	4 cm then 5 cm
	D	8.25 ± 5.92	61.7	45.6	95.6
	D x 2	6.97 ± 5.26	68.2	54.1	92.2
	D x 4	5.17 ± 4.07	76.3	67.3	88.0
	D x 6	4.17 ± 3.38	80.8	77.5	84.4
	D x 8	3.85 ± 3.21	82.3	84.8	80.1
	D x 10	4.15 ± 3.38	81.3	89.3	74.6
	D x 12	4.88 ± 3.87	78.6	92.0	68.6
	D x 14	5.83 ± 4.53	75.0	93.6	62.6

Post-op lesion Computed lesion Simulation with nominal parameters Simulation with personalized parameters Portal vein Tumor Computed lesion

Table 5 .

 5 1: Ranges of parameters values explored in the sensitivity analysis. lesion obtained with nominal values. Thoses values are reported in Table5.2. As it is not reasonable to try to estimate all of these four parameters at once, we decided to estimate only two of them: we chose dt for its effect on the volume, as the nominal value of kf [ONeill 2011] seemed accurate and ε. H was choosen large enough to maintain a constant temperature of 37 • C in the CT-visible vessels. We fixed the other parameters to nominal values for the personalization on patient data.

			Total Sobol indices
	Parameters volume	point-to-mesh error
	dt	0.58	0.37
	H	0.15	0.35
	kf	0.16	0.16
	ε	0.43	0.37

Table 5 .

 5

2: Total Sobol indices on the volume of the lesion and the average point-tomesh error with respect to lesion obtained with nominal parameters.

Table 5 .

 5

	Estimated Parameters

3: Evaluation of the parameter estimation framework on patient data: 7 tumors from 3 patients are considered.

Temperature Profile during heating phase Temperature Profile during cooling phase Power Profile during heating phase

	Pre-clinical Data	Pre-processing	RFA simulation
	CT of the probe	Probe Segmentation	Estimation of the Probe Position	RFA computation:
	Anatomical MRI with RFA probe	Pre-op CT-MRI registration	Estimation of Liver Anatomy from CT	Temperature evolution +
	Pre-op CT Images	Segmentations		Delivered Power + Cell State
	Phase-Contrast MRI Invasive catheter Probe Thermistors	Blood Flow Measurements Pressure Measurements	Computation of 3D blood Flow	Evolution during the RFA simulation + -Computed Temperature Profile	Parameter Estimation
	Probe Measurements			+ -
				Computed Power
	Post-op CT or	Segmentation of post-		Profile
	MRI Images	op necrotic regions		Computed
				Necrotic regions
	Post-op data Intra-op data Pre-op data	Post-op to Pre-op registration		+	-

Table 6 .

 6 1: Nominal parameter values.

	Notation	Parameter Name	Nominal
	T k ( • C)	parameter of cell death model	40.5
	β (s -1 )	damage rate coefficient	3

Table 7 .

 7 2: Evaluation of the smooth vessel trees: Dice, Positive Predictive Value (PPV) and sensitivity between the original and the smoothed vessel images are computed for Pig 4.

	Mesh	DICE (%)	PPV (%)	Sensitivity (%)
	arteries	93.1	93.9	92.3
	portal vein	91.3	85.3	98.3
	hepatic veins	84.9	93.8	77.5
	DICE, Positive Predictive Value (PPV) and sensitivity scores have already been
	defined in Chapter 4 in Section 4.4.1, but we will repeat the definition here for the
	sake of clarity.

Table 7 . 4 :

 74 Comparisons between personalized and nominal values of estimated parameters.

	Notation	Nominal	Personalized Values	Automatic
		Value	on tumor 1	on tumor 4-3	optimization
	δ	3.3×10 -3	1×10 -4	1×10 -4	No
	c U t dt	3.6 ×10 3 0.512	3.6 ×10 1 0.614	3.6 ×10 1 0.512	Yes Yes

Table 7 .

 7 7: Quantitative evaluation of the probe position on Pig 3-2: on average, point-to-mesh errors decrease significantly and DICE scores increase by a factor 2.

	Pig 3-2	Power error	Temperature	point-to-mesh error in	Dice	Sensitivity	PPV
	Probe in	(W)	error ( • C)	mm (mean, sd, min,	(%)	(%)	(%)
				max, median, P 95 )

Table 7 .

 7 8: Evaluation of the effect of the cell death model -thermal Dose at 240 on the point-to-mesh errors and DICE scores between the simulated lesion and ground truth on Pig 2-2.

	Pig 2-2	point-to-mesh error in mm	Dice	Sensitivity	PPV
		(mean, sd, min, max, median,	(%)	(%)	(%)
		P 95 )

Table 7 .

 7 9: Quantitative evaluation of model prediction: Effect of the cell death model on the shape of the simulated lesion.

	Pig 2-2	sphericity	volume	max diameter
			(mm 3 )	(mm)
	(GT)	0.45	12853	50.09
	(GT-smooth)	0.71	5356	30.43
	Thermal Dose
	(nominal-0.5)	0.96	12519	34.19
	(nominal-0.1)	0.93	17019	37.39
	(perso)	0.94	14641	35.37
	Cell Death
	(nominal-0.5)	0.90	5406	23.96
	(nominal-0.1)	0.90	7193	26.32
	(perso)	0.90	5933	24.78

Table 7 .

 7 10: Quantitative evaluation of model prediction. The mean errors between simulated and measured values are computed for each ablation. For the eight ablations which are performed through several cooling and heating periods, the different error values for each period are reported. For each tumor, several computations are performed, with nominal parameters and two different time-step (nominal-dt), with personalized parameters, and with advection. The parameter estimation has been performed on Pig 1 and tumor 4-3 and gave a set of personalized parameters, which are used in every case (perso). When it is not specify, the time-step is dt = 0.1 s. For computation with nominal parameters (respectively with personalized parameters), on average, mean errors of 25.9 W (resp. 20.4 W), 5.1 • C (resp. 7.8 • C) and 5.32 mm (resp. 5.99 mm) of mean of the point-to-mesh error are achieved. This table presents results for Pig 1 and Pig 4.

	Pig	Power error	Temperature	point-to-mesh error	Dice	Sensitivity	PPV
		(W)	error ( • C)	in mm	(%)	(%)	(%)
				(mean, sd, min,
				max, median, P 95 )
	1
	(nom-0.5)	21.9	8.3	5.14, 4.77, 0.14,
				21.38, 3.51, 16.13

Table 7 .

 7 11: Quantitative evaluation of model prediction. This table presents results for Pig 2.

	Pig	Power error	Temperature	point-to-mesh error	Dice	Sensitivity	PPV
		(W)	error ( • C)	in mm	(%)	(%)	(%)
				(mean, sd, min,
				max, median, P 95 )
	2-1
	(nom-0.5)	7.4	2.6	2.93, 2.66, 0.18,	59.8	44.1	92.8
				15.83, 2.05, 8.88
	(nom-0.1)	16.0	1.5	2.46, 2.28, 0.09,	66.9	54.3	87.0
				14.65, 1.64, 7.06
	(perso)	10.6	1.3	2.81, 2.55, 0.12,	62.0	46.8	92.0
				15.25, 1.97, 8.42
	2-2
	(nom-0.5)	47.5	3.1	8.52, 5.03, 0.11,	13.8	9.8	23.3
				23.30, 8.37, 16.59
	(nom-0.1)	62.3	5.0	8.03, 4.85, 0.10,	16.8	13.1	23.5
				22.28, 7.83, 15.58
	(perso)	50.0	4.5	8.44, 5.00, 0.10,	14.3	10.4	22.6
				22.99, 8.27, 16.35

Table 7 .

 7 12: Quantitative evaluation of model prediction. The tumor 3-1 was ablated twice with two consecutive heating periods, thus the simulated necrosis cannot be compared with any ground truth. This table presents results for Pig 3.

	Pig	Power error	Temperature	point-to-mesh error	Dice	Sensitivity	PPV
		(W)	error ( • C)	in mm	(%)	(%)	(%)
				(mean, sd, min,
				max, median, P 95 )

Table 7 .

 7 13: Quantitative evaluation of model prediction. This table presents results for Pig 5.

	Pig	Power error	Temperature	point-to-mesh error	Dice	Sensitivity	PPV
		(W)	error ( • C)	in mm	(%)	(%)	(%)
				(mean, sd, min,
				max, median, P 95 )
	5-1
	(nom-0.5)	9.6, 23.0, 25.2,
		28.8

Table 7 .

 7 14: Quantitative evaluation of model prediction on the shape of the computed lesion. Results for Pig 1, Pig 2 and Pig 3 are reported here.

	Pig	sphericity	volume	max diameter
			(mm 3 )	(mm)
	1 (GT)	0.71	10228	43.46
	1 (nominal -0.5)	0.83	4137	27.48
	1 (nominal -0.1)	0.84	5450	29.27
	1 (perso)	0.83	4423	28.05
	1 (advection)	0.82	8568	35.86
	2-1 (GT)	0.77	13337	42.67
	2-1 (nominal-0.5)	0.90	6344	24.41
	2-1 (nominal-0.1)	0.90	8316	26.65
	2-1 (perso)	0.90	6781	24.88
	2-2 (GT)	0.45	12853	50.09
	2-2 (nominal-0.5)	0.90	5406	23.96
	2-2 (nominal-0.1)	0.90	7193	26.32
	2-2 (perso)	0.90	5933	24.78
	3-1 (GT)	N.A	N.A	N.A
	3-1 (nominal-0.5)	0.84	13662	33.76
	3-1 (nominal-0.1)	0.83	17616	37.22
	3-1 (perso)	0.84	16004	36.18
	3-1 (advection)	0.83	17593	38.69
	3-2 (GT)	0.61	18105	57.80
	3-2 (nominal-0.5)	0.87	46048	48.96
	3-2 (nominal-0.1)	0.86	55386	52.44
	3-2 (perso)	0.86	53008	51.79
	3-3 (GT)	0.69	21722	45.18
	3-3 (nominal-0.5)	0.89	13059	32.89
	3-3 (nominal-0.1)	0.89	17370	35.64
	3-3 (perso)	0.89	15781	34.66
	3-3 (advection)	0.87	19934	39.53

Table 7 .

 7 15: Quantitative evaluation of model prediction on the shape of the computed lesion. Results for Pig 4 and Pig 5 are reported here. to enable parameter estimation, to test the generalization of a biophysical model in terms of tissue properties, but also for the expected clinical applications. For interventional training, interactive therapy planning or procedure guidance, faster than real-time computations would be ideal. We have not reached this point yet, but with our LBM implementation, near real-time computation can be achieved on a standard machine.

	Pig	sphericity	volume	max diameter
			(mm 3 )	(mm)
	4-1 (GT)	0.66	1849	23.94
	4-1 (nominal-0.5)	0.89	10956	31.70
	4-1 (nominal-0.1)	0.89	14436	34.73
	4-1 (perso)	0.89	12574	33.03
	4-2 (GT)	0.84	13066	35.86
	4-2 (nominal -0.5)	0.88	4001	22.98
	4-2 (nominal -0.1)	0.89	5621	24.94
	4-2 (perso)	0.88	4755	23.77
	4-3 (GT)	0.86	11322	35.19
	4-3 (nominal-0.5)	0.86	8466	30.33
	4-3 (nominal-0.1)	0.85	10778	32.65
	4-3 (perso)	0.85	9441	31.50
	5-1 (GT)	0.80	15849	42.59
	5-1 (nominal-0.5)	0.88	12667	32.45
	5-1 (nominal-0.1)	0.88	16040	35.19
	5-1 (perso)	0.88	14486	33.97
	5-1 (advection)	0.88	16040	35.19
	5-2 (GT)	0.71	19129	51.72
	5-2 (nominal-0.5)	0.80	6339	28.79
	5-2 (nominal-0.1)	0.79	8562	31.62
	5-2 (perso)	0.79	7552	30.43
	5-2 (advection)	0.79	8398	31.21
	5-3 (GT)	0.69	7432	35.89
	5-3 (nominal-0.5)	0.85	5237	24.70
	5-3 (nominal-0.1)	0.85	6939	26.76
	5-3 (perso)	0.85	6001	25.98
	5-3 (advection)	0.85	6799	26.76

•

  Chloé Audigier, Parameter Estimation For Personalization of Liver Tumor Radiofrequency Ablation. Siemens Corporate Research, December, 4th, 2014, Princeton, USA.

	Awards
	• Outstanding Paper Award at the 6th International Workshop on Abdom-
	inal Imaging: Computational and Clinical Applications, 2014.

  This equation states that the total rate of change of the distribution function is equal to the rate of the collision. The total rate of change can be expressed as:

	df =		∂f ∂r	dr +	∂f ∂c	dc +	∂f ∂t	dt
	Dividing by dt:
	df dt	=	∂f ∂r	dr dt	+	∂f ∂c	dc dt	+	∂f ∂t
	And:
	df dt	=	∂f ∂r	c +	∂f ∂c	a +	∂f ∂t

  ∇∇ : t e i .e i (f i )) + 0(∆t) 3where the operator : represents tensor contraction. Introducing scaling, keeping term up to ε 2 we get:f i (p + ce i ∆t, t + ∆t) = f i (p, t) + ε ∂f i ∂t 1 ∆t + ε 2 ∂f i ∂t 2 ∆t + εe i .∇f i c∆t (ce i .∇f i ) + ε 2 c 2 ∇∇ : f i t ei .e i ) + 0(∆t) 3

				t) +	∂f i ∂t	∆t + ∆tce i .∇f i +	∆t 2 2	(	∂ 2 f i ∂t 2 + 2	∂ ∂t	(ce i .∇f i )
			+c 2 + ∆t 2 2 (ε 2 ∂ 2 f i ∂t 2 1 + 2ε 2 ∂ ∂t 1
		Substituting this equation in equation A.6:
	ε	∂f i ∂t 1	∆t + ε 2 ∂f i ∂t 2	∆t + εe i .∇f i c∆t +	∆t 2 2	(ε 2 ∂ 2 f i ∂t 2 1	+ 2ε 2 ∂ ∂t 1	(ce

i .∇f i ) + ε 2 c 2 ∇∇ : f i t e i .e i ) + 0(∆t) 3 = w[f eq i (p, t) -f i (p, t)]

  t) + ∂f i ∂t ∆t + ∆tce i .∇f i .∇f i ) + c 2 ∇∇ : t e i .e i (f i )) + 0(∆t) 3where the operator : represents tensor contraction.Introducing scaling, keeping term up to ε 2 we get:f i (p + ce i ∆t, t + ∆t) = f i (p, t) + ε ∂f i ∂t 1 ∆t + ε 2 ∂f i ∂t 2 ∆t + εe i .∇f i c∆t (ce i .∇f i ) + ε 2 c 2 ∇∇ : f i t e i .e i ) + 0(∆t) 3 (ce i .∇f i ) + ε 2 c 2 ∇∇ : f i t e i .e i ) + 0(∆t) 3 = w[f eq i (p, t) -f i (p, t)] + ∆tω i (εS 1 + ε 2 S 2 )

			+ (ce i + ∆t 2 2 ( ∂ 2 f i ∂t 2 + 2 ∂ ∂t ∆t 2 2 (ε 2 ∂ 2 f i ∂t 2 1 + 2ε 2 ∂ ∂t 1
		Substituting this equation in equation A.8:
	ε	∂f i ∂t 1	∆t + ε 2 ∂f i ∂t 2	∆t + εe i .∇f i c∆t +	∆t 2 2	(ε 2 ∂ 2 f i ∂t 2 1	+ 2ε 2 ∂ ∂t 1

  So by adding the first equation times ε and the second one times ε 2 , we get:There is an error term:ε 2 (τ -∆t 2 ) ∂ ∂t 1 [∇.(uT + S 1 )] . The relaxation parameter can be related to the diffusion coefficient as:We use Neumann Boundary condition, with an assumption of zero flux through the boundary. Thus, if r is on the boundary and r + ce α ∆t is outside the domain, then: + ∆t) = 2∆[ fα (r, t + ∆t)] + (1 -2∆)[ fα (r + ce ᾱ∆t, t + ∆t)]

	And:								N
									f 0 i	t e i .e i = T c 2 s I
									i=0
	So we get:
		∂T ∂t 2	+ (	∆t 2	-τ )(	∂S 1 ∂t 1	+	∂ ∂t 1	[∇.(uT )] + c 2 c 2 s ∇ 2 T ) = S 2
	Finally:
		∂T ∂t 2	+ c 2 c 2 s (	∆t 2	-τ )∇ 2 T + (	∆t 2	-τ )	∂ ∂t 1	[∇.(uT ) + S 1 ] = S 2
	∂T ∂t	+ u.∇T = c 2 c 2 s (τ -	∆t 2	)∇ 2 T + ε(τ -	∆t 2	)	∂ ∂t 1	[∇.(uT ) + εS 1 ] + R
									d = (τ -	∆t 2	)c 2 s c 2
	A.5 Boundary Condition
	for ∆ ≥ 1 2
		f ᾱ(r, t + ∆t) =	1 2∆	[ fα (r, t + ∆t)] +	2∆ -1 2∆	[ fᾱ (r, t + ∆t)]
	for ∆ < 1 2
	f ᾱ(r, t
		∂f 0 i ∂t 2	+ (	∆t 2	-τ )(	∂S 1 ∂t 1	+	∂ ∂t1	(

N i=0 ce i .∇f 0 i ) + c 2 ∇∇ : N i=0 f i t e i .e i )) = S 2

www.cancer.org/the-history-of-cancer-pdf Retrieved

2015-07-15.

Radiofrequency ablation. (2015, July 13). In Wikipedia, The Free Encyclopedia. Retrieved 16:44, August 14, 2015, from https://en.wikipedia.org/w/index.php?title=Radiofrequency_ ablation&oldid=671315922

High-intensity focused ultrasound. (2015, May 4). In Wikipedia, The Free Encyclopedia. Retrieved 14:08, August 5, 2015, from https://en.wikipedia.org/w/index.php?title= High-intensity_focused_ultrasound&oldid=660711436

 7 Microwave ablation.(2014, May 10).In Wikipedia, The Free Encyclopedia. Retrieved 14:10, August 5, 2015, from https://en.wikipedia.org/w/index.php?title=Microwave_ ablation&oldid=607944914

http://dakota.sandia.gov -multilevel framework for sensitivity analysis

http://www.paraview.org/

The Vascular Modeling Toolkit website, www.vmtk.org.

www.cgal.org.

www.sofa-framework.org.

http://dakota.sandia.gov -multilevel framework for sensitivity analysis.

Chapter 7

Pre-Clinical Study: Results and Discussion

Ground truth

computational model of RFA and to consider potential safety issue of the proposed application. Even if promising results are achieved with the use of patient-specific parameters, the impact of possible biases in the post-to the pre-operative image registration like the impact of the average simulation of the probe need to be investigated as well as the sensitivity of the results with respect to segmentation and registration.

Chapter 8

Summary and Perspectives

Essentially, all models are wrong, but some are useful. Each chapter of this manuscript contains conclusions and perspectives on the performed work. We now conclude on the main achievements of this PhD work. Then, we suggest possible improvements related to this work.

Main Contributions

Modeling of Radiofrequency Ablation with LBM

We have presented an integrated framework for fast modeling of Radiofrequency Ablation (RFA), which computes and predicts RFA interventional outcome in terms of necrotic zone, delivered power and temperature distribution throughout the liver. Our detailed computational model of the biophysical mechanisms (heat transfer, cellular necrosis, hepatic blood flow) involved in RFA of abdominal tumors is based on patient images. We focused on modeling heat propagation and cellular necrosis based on a patient image while considering the heat sink effect of blood vessels and porous circulation in the liver. It is implemented on Graphics Processing Units (GPU) using the Lattice Boltzmann Method to take advantage of the computation capacity of the Graphics Card and thus to reach near real-time computation. Fast

List of Publications

The presented work led to several published and submitted publications as first author.

Appendix A

The Lattice Boltzmann Method

A.1 The Boltzmann Equation

A.1.1 Boltzmann Transport Equation where 0 and I are respectively the zero matrix and identity matrix. The value of c s vary from one configuration to another (depends on the number of lattices). The last equation can also be written as :

where δ α,β is the Kronecker's delta.

The last isotropic property is:

ω i e i,α e i,β e i,γ = 0

The first momentum equation is:

f eq i By taking the sum, we get:

The second momentum equation gives:

f eq i ce i

So:

By expanding the scalar product: ω i e i,α e i,β

According to the last isotropic properties, we get: