Keywords: photodynamic therapy, glioblastoma, dosimetry, Monte-Carlo simulation

Photodynamic therapies of high-grade gliomas: from theory to clinical perspectives

Gliomas are the most common primary brain tumors in adults. Among them, glioblastoma (GBM) represents the most frequent primary brain tumor and have the most dismal prognosis. Its annual incidence is about 3 to 5 cases for 100,000 persons (about 3000 news cases each year in France). Median survival varies between 14 to 15 months according to the extent of tumor resection.

The standard of care includes surgery and is followed by radiation therapy and chemotherapy. Maximal resection is expected to delay recurrence. Despite of using intraoperative photodynamic diagnosis, or fluorescence guided resection (FGR), which improves the extent of resection, relapse still occurs in these resection margins in 85% of cases.

Alternatives therapies have to be developed to enhance patients' overall survival. In this context, Photodynamic Therapy (PDT) seems relevant. PDT is based on the synergy of three parameters: a photosensitizing molecule, the photosensitizer (PS) that concentrates preferentially into the tumor cells, laser light and oxygen. Laser light induces a reaction between the PS and the oxygen of the cell. This reaction produces highly cytotoxic molecules (including singlet oxygen) and leads to death of tumor cells. Two treatment modalities are investigated: interstitial PDT (iPDT) or intraoperative PDT.

The main goal of this thesis is to provide technological tools to develop the PDT for GBM treatment. Thus, the two treatment modalities have been investigated.

When tumor resection is non-achievable (approximately 20% to 30% of cases), iPDT may be preferred. This modality aims to insert optical fibers directly into the target to illuminate tumor tissues. Thus, simulation of light propagation in brain tissues is required to plan the location of optical fibers. Considered as reference method, a Monte-Carlo model accelerated by graphics processing unit was developed. This model computes the light propagation emitted by a cylindrical diffusor inside heterogeneous media. Accuracy of the model was evaluated with experimental measurements. The acceleration provided by the parallelization allows its use in clinical routine.

The iPDT has to be planned using a Treatment Planning System (TPS). A proof of concept of a TPS dedicated to the stereotactic iPDT treatment of GBM was developed. This software provides basic tools to plan the stereotactic insertion of cylindrical diffusors in patient's brain and to compute the associated dosimetry. The stereotactic registration and the dosimetry computation's accuracy were evaluated with specific methodologies.

When tumor resection is achievable, the intraoperative PDT may be applied early after the FGR. It takes advantage of the presence of the PS (protoporphyrin IX) used for FGR purpose and that is already concentrates into the tumor cells. Thus, the proposed treatment strategy fits into the current standard of care. A medical device was designed into fit to the resection cavity and illuminate homogeneously the cavity's margins. The device is constituted of two parts: a trocar coupled to an inflatable balloon and a fiber guide developed in the ONCO-THAI laboratory allowing to insert the light source. Specific methodologies were developed to calibrate and assess the device in terms of mechanical properties and dosimetry.

The calibration process leaded to a transfer function that provides fast, robust and easy treatment duration prescription to induce a PDT response in cavity margins. Furthermore, a comprehensive experimental design has been worked out prior to the clinical trial that evaluate the safety of the procedure.

Thérapies photodynamiques appliquées aux gliomes de haut grade : de la théorie à la réalité clinique

Les gliomes sont les tumeurs cérébrales primaires les plus communes chez l'adulte.

Parmi eux, le glioblastome (GBM) représente la tumeur cérébrale la plus fréquente avec le pronostic le plus sombre. Son incidence annuelle est d'environ 3 à 5 cas pour 100 000 personnes (environ 3000 nouvelles chaque année en France). La survie médiane varie entre 14 et 15 mois selon la qualité de la résection tumorale.

Le standard de soins inclut une résection chirurgicale suivie d'une radiothérapie et d'une chimiothérapie. Une résection maximale est souhaitée afin de diminuer les risques de récidive. Bien que l'utilisation de la technique de diagnostic photodynamique peropératoire, appelée résection fluoroguidée (FGR), améliore la qualité de résection, une récidive survient dans ces berges de la cavité opératoire dans 85% des cas.

Des thérapies alternatives doivent être développées pour améliorer la survie globale des patients. Dans ce contexte, la thérapie photodynamique (PDT) semble pertinente. La PDT est basée sur la synergie de trois paramètres : une molécule, le photosensibilisateur (PS) qui se concentre préférentiellement dans les cellules tumorales, la lumière laser et l'oxygène. La lumière laser induit une réaction entre le PS et l'oxygène de la cellule. Cette réaction produit des molécules cytotoxiques (dont l'oxygène singulet) et conduit à la mort de cellules tumorales. Deux modalités de traitement sont étudiées : la PDT interstitielle (iPDT) ou la PDT peropératoire.

L'objectif principal de cette thèse est de fournir des outils technologiques afin développer la PDT pour le traitement du GBM. Ainsi, les deux modalités de traitement ont été étudiées.

Lorsque la résection n'est pas réalisable (environ 20% à 30% des cas), l'iPDT peut être privilégiée. Cette modalité vise à insérer des fibres optiques dans la cible thérapeutique pour éclairer les tissus tumoraux. Ainsi, la simulation de la propagation de la lumière dans les tissus est nécessaire pour planifier la localisation des fibres optiques. Considéré comme méthode de référence, un modèle Monte-Carlo accéléré par processeurs graphiques a été développé. Ce modèle calcule la propagation de la lumière émise par un diffuseur cylindrique dans des milieux hétérogènes. La précision du modèle a été évaluée avec des mesures expérimentales. L'accélération fournie par la parallélisation permet son utilisation dans la routine clinique. L'iPDT doit être planifiée à l'aide d'un système de planification de traitement (TPS). Une preuve de concept d'un TPS dédié au traitement stéréotaxique iPDT du GBM a été développée. Ce logiciel fournit des outils de base pour planifier l'insertion stéréotaxique de diffuseurs cylindriques et calculer la dosimétrie associée. Le recalage stéréotaxique et la précision du calcul dosimétrique ont été évalués avec des méthodologies spécifiques.

Lorsque la résection est réalisable, la PDT peropératoire peut être appliquée juste après la FGR. Celle-ci profite de la présence du PS (la protoporphyrine IX) utilisé pour la FGR et qui s'est déjà concentrée dans les cellules tumorales. Ainsi, la stratégie de traitement proposée peut s'inclure facilement au standard de soin. Un dispositif médical a été conçu pour s'adapter à la cavité et éclairer de façon homogène les berges de la cavité opératoire. Le dispositif est constitué de deux parties : un trocart couplé à un ballon gonflable et un guide de fibre optique développé au sein du laboratoire ONCO-THAI permettant d'insérer la source lumineuse. Des méthodologies spécifiques ont été développées pour étalonner et évaluer l'appareil en termes de contraintes mécaniques et dosimétrique. L'étalonnage a permis la création d'une fonction de transfert permettant une prescription de durée de traitement rapide, robuste et facile. De plus, de nombreux tests ont été réalisés en amont de l'essai clinique qui évalue la sécurité de la procédure.

Mots clés : thérapie photodynamique, glioblastome, dosimétrie, simulation Monte-Carlo, calcul parallèle par processeur graphique, dispositif médical, étude clinique

General introduction

Context

Gliomas are the most common primary brain tumors in adults. Among them, glioblastoma (GBM) has the most dismal prognosis. It represents the most frequent primary brain tumor in adults [1][2][3] with an annual incidence of about 3 to 5 cases for 100,000 persons [4][5][6]. Even if no environmental risks factor has been identified except in radiation therapies [7,8], the population growth and aging, the diagnostic improvement and the dismal prognosis of this tumor are expected to challenge the research community to find relevant answers to this medical emergency [9]. The first-line treatment is usually surgery, either to confirm the diagnosis with a biopsy or to remove as much of the tumor as possible. Once the surgery or biopsy has been realized, adjuvant treatments, such as radio and/or chemotherapy, are administered according to the grade and the genetic profile of the tumor. Unfortunately, this standard of care has a modest result regarding the patient's survival. Thus, alternative therapies using different ways of treatment are expected to bring effective solutions to GBM treatment.

Among the recent alternatives, the photodynamic therapy (PDT) may be a promising treatment in oncology and particularly in the treatment of brain cancers, including GBM. PDT is a non-thermal energy-based therapy relying on light exposure after photosensitization of tumor cells. A PDT effect is the result of the synergy of three elements: a photosensitizing molecule, the energy delivered through laser light exposure in order to excite the molecule, and the oxygen. The combination of these three components leads to the formation of reactive oxygen species, including singlet oxygen, which are cytotoxic molecules that damage the tumor cells.

PDT was firstly used in the Antiquity to treat skin disease [10] and was reconsidered in 1972 [11]. The first clinical trial using PDT on GBM was realized in 1980 by Perria et al. [12]. The advent of this therapy is largely related to the improvement of laser systems, particularly in medicine. Currently, two clinical PDT modalities for neurosurgery are investigated. When the surgical resection cannot be achieved because of surgical constraints, the PDT is delivered interstitially. This treatment modality plans to insert optical fibers under stereotactic conditions and to vehicle light directly into the tumor. When the patient can undergo to surgery, PDT is delivered intraoperatively, i.e., early after the resection. This lighting modality aims at illuminating the cavity created by the surgical resection in order to treat the remaining tumor tissue that cannot be removed.

Although there are few studies demonstrating proofs of concepts, this therapy is not adapted to current clinical therapeutic modalities for several issues yet. No consensus has been defined on the use of a photosensitizer and there is a lack of reliable and reproducible therapy delivery systems (medical devices and treatment planning software). All these disparities lead to the absence of randomized, controlled multicenter clinical trials on the treatment of GBM by PDT.

Contributions

This work has been achieved in the ONCO-ThAI laboratory (Inserm, University of Lille, Hospital of Lille), which develops therapies based on laser light, mostly in oncology.

An interdisciplinary research is led in order to reach the clinic evaluation of the technologies developed in the laboratory.

The main objective of this thesis was to develop technological tools adapted to clinical routine, both software and hardware, dedicated to the treatment of GBM by PDT. Several axis of research have been investigated in which the dosimetry was the guideline: planning the PDT in both lighting modalities (interstitial and intraoperative) and designing the device to deliver the treatment.

The main contributions of this thesis can be summarized as follows:

-A GPU-accelerated Monte-Carlo method is proposed to compute the light distribution emitted from cylindrical diffusers in heterogeneous biological tissues.

-A proof of concept of a treatment planning system dedicated to stereotactic interstitial photodynamic treatment is introduced. The workflow of this software is similar to the one used in neuro-radiosurgery. It proposed to plan the stereotactic insertion of several cylindrical diffusers and compute the light propagation during the treatment procedure.

-A new medical device dedicated to intraoperative photodynamic therapy has been designed and characterized using ex vivo experimentations. Results have been compared to Monte-Carlo simulations.

-A transfer function has been created in order to easily obtain the treatment duration according to the size of the cavity to obtain a therapeutic effect at 5 mm inside margins resection.

-A clinical trial has been set up to evaluate the safety of the procedure using this medical device. This phase one clinical trial is currently under evaluation.

Roadmap

This dissertation is organized as follows.

Chapter I introduces laser technologies and their application in medicine including the PDT. A focus is made on the photosensitizer and the biological mechanisms involved. Then, a description of the glioblastoma, its standard of care and alternative treatments are presented.

Among them, previous photodynamic treatments of GBM, including interstitial and intraoperative modalities treatment, are listed.

In order to plan PDT treatments, a dosimetry modelling is needed. Since light is one of the three major components of PDT, Chapter II details the physics interactions of light propagation in biological tissues used in the predictive dosimetry. Among all dosimetric features, models of light delivery computing are presented. Also, some monitoring methods used to evaluate the PDT dosimetry are described.

The Monte-Carlo method used to model the light propagation in tissue is still the reference model. Chapter III presents this model and techniques to reduce the computing time.

Among them, a Monte-Carlo model based on acceleration by GPU parallel computing is proposed. Its implementation and the evaluation methodology are presented. Also, optimization of algorithm's performances is described.

A proof of concept of a treatment planning system dedicated to stereotactic interstitial photodynamic treatment is detailed in chapter IV. Methodologies to validate the stereotactic registration and evaluate the dosimetry part are described. Future improvements including the integration of the GPU Monte-Carlo model presented in chapter III are discussed.

Finally, Chapter V is dedicated to the intraoperative PDT treatment. The design of a new lighting applicator is described. Methodologies to characterize the device is detailed, including its dosimetry, a mechanical assessment and other safety tests. Thus, the method proposed enables to obtain the treatment duration according to the cavity size in order to obtain a therapeutic effect in resection margins. The clinical investigation of the device within the clinical trial INDYGO is described. Further improvements and future technologies to substitute this medical device are discussed. This chapter presents the concept of photodynamic therapy (PDT) dedicated to the treatment of High-Grade Gliomas (HGG). An overview on laser is firstly given which includes a brief history of laser technology, the physics principles, the main technological components and its current applications in medicine. Then, photosensitizers are presented with a specific focus on the precursor, 5-aminolevulinic acid . Direct and indirect processes of tumor cells destruction resulting from PDT are highlighted. In the context of this study, primary brain tumors, especially glioblastoma (GBM) and their standard of care are described before a systematic review of PDT technics dedicated to treat HGG.

I. Photodynamic therapy 1. LASER: History

"This property of elementary processes as expressed by equation (12) makes a quantum theory of radiation almost unavoidable. The weakness of the theory lies, on the one hand, in its not bringing us closer to a union with the wave theory, and, on the other hand, that it leaves the time and direction of the elementary processes to chance; in spite of this, I have full confidence in the trustworthiness of this approach." [13] This quote comes from a paper published by Albert [START_REF] Einstein | On the quantum mechanics of radiation[END_REF] In this paper, the scientist confirms his interest in the Planck's theory according to which light radiation is composed of several "quanta" (namely photons) that can be absorbed or emitted discontinuously. Furthermore, he carries on the corpuscular theory of light and its duality with the wave theory and describes for the first time the stimulated emission of radiation process [14]. Later, this physical principle will be used in the MASER technology (Microwave Amplification by Stimulated Emission of Radiation), the ancestor of the well-known LASER (Light Amplification by Stimulated Emission of Radiation) described by Gould in 1959 [15] and firstly designed by Maiman in 1960 [16, 17].

a. Physics principles

The laser is based on two main physic processes: the stimulated emission of radiation and the population inversion.

The stimulated emission of radiation occurs when a photon with a particular energy encounters an excited atom (see figure 1.1). If the energy level hv of the incident photon is equal to the difference between levels of the electrons E1 and E0, the incident photon then stimulates the electron return to the ground state E0. This relaxation leads to the creation of a second photon bearing the same properties of the incident photon. Regarding laser, light beam of 1J at a wavelength of 630nm requires about 3.10 19 photons using stimulated emission process. If the two photons resulting from this stimulated emission process do not meet other excited atoms, then, they have a strong probability to be absorbed. It appears thus necessary to have a majority of excited atoms to produce a large number of photons and because atoms are mostly naturally at the ground state, another energy source is expected to excite atoms.

A second issue occurs if we consider only two states atoms constituting the irradiated media. To create enough photons from stimulated emission process, photons must have a higher probability to meet an excited atom than to meet a stable one. However, if a photon is emitted in a media composed of 50% of excited atoms at an energy level E1 and 50% of ground state atoms E0, statistically, only absorption and spontaneous emission processes is obtained (independently from the initial conditions of the atoms population). Thus, a third state must be considered to perform a "population inversion" (see figure 1.2). This term defines the process during which a media owns more excited atoms than natural ground state atoms. This inversion can occur only with at least a triple states system.

Let us consider a media with triple states atoms: at the Boltzmann equilibrium, there are more atoms at the ground state E0 than excited state E1, and there are more atoms at the excited state E1 than excited state E2. However, if an external source of energy Eext equal to E2 -E0 is applied to this system, several E0 atoms are promoted to the E2 energy level. The number of E2 atoms tends towards the same number of E0 atoms but cannot exceed this number. Nonetheless, the population of E2 atoms can be larger than the population of E1 atoms. In this case, the inversion population between E2 and E1 is reached. An incident photon owns therefore more chances to encounter an excited atom and to provoke a stimulated emission corresponding to a relaxation of E2 -E1. The stimulated emission of two excited atoms by photons from a first stimulated emission induces four photons and consequently leads to a cascade of photons. If the excited atoms proportion is maintained higher than the proportion of ground state atoms and if the photons created by stimulated emission remain with a sufficient number for them to induce stimulated emission processes, a laser light is created. Mainly due to the conservation of photon's properties during the stimulated emission process, each photon composing the light beam owns the same features.

Thus, laser light is:

-Monochromatic: each photon has the same wavelength. However, due to several mechanisms of spectral emission broadening, a laser light cannot be a pure monochromatic light beam.

-Directional: the light beam produced by a laser has a very small divergence in space. This strong directionality allows to obtain a highly focal beam.

-Temporally coherent: because all photons produced by a laser are almost with the same wavelength, each photon's phase is synchronized and is the same for all of them.

-Spatially coherent: additionally, in a perpendicular plan to the light propagation, all photons have the same electromagnetic field (value and direction).

Using these particular properties, laser light has quickly found several applications in many different fields.

b. Technological setting

To produce stimulated emission and population inversion, all lasers are composed of three different parts: active gain medium, pumping and resonating cavity.

Active gain medium is the core of a laser, where the stimulated emissions happen. It is composed of atoms with determined energy levels that define the wavelength of the laser light. Changing the composition of this part can modify the wavelength of the laser light.

Several active gain media are employed and allow the classification of the laser type according to the material composition:

-Solid state: crystal (ruby, Nd-YAG) and optical fiber (silicate or phosphate glasses) -Semi-conductor material -Liquids: dyes -Gases (CO², He-Ne, argon) Amplifier or laser pumping is dedicated to the realization of the inversion population.

It transfers an external energy to the active gain medium in order to excite atoms and to produce in return stimulated emission of radiation. Several types of pumping have been developed including:

-Optical pumping using flash lamp, arc lamp or the use of an additional laser -Electrical pumping: application of an external electric current or electron beam, glow discharge -Other types such as thermic pumping (adiabatic cooling), chemical reaction and even using nuclear reaction products.

Resonating cavity is composed of two face-to-face mirrors and one is partially transparent. The proportion of transparency must be highly controlled. If too many photons leave the resonator, the percentage of stimulated emission decreases and tends toward zero.

Furthermore, the distance between those two mirrors represents one of the main filters to obtain an output light beam with a narrow waist of the spectral emission peak. The distance is computed in order to attenuate most of non-desired wavelengths using resonance phenomenon.

c. Applications

Decades ago, laser device was highly unstable, expensive, cumbersome, with a very bad efficiency and was mostly dedicated to research applications. With the advent and the development of the semi-conductor technology, laser technology has been strongly democratized in many applications, including medicine. Nowadays, laser systems in medicine are more efficient, portable and can be highly powerful (see figure 1.3). Most of lasers currently manufactured are based on laser diode technology. Its three parts (active gain medium, pumping and resonating cavity) are contained in an electronic component of only few millimeters. Output power can easily reach several watts. Additionally, with the improvement in the electronic doping and recent semi-conductor compositions, a large range of wavelength is possible (from 405nm with Indium-Gallium-Nitride compound to 980nm with Indium-Gallium-Arsenide compound). All these improvements have led to a larger use of the laser in the medical field.

Photosensitizer

The use of molecules for which the cytotoxic effect can be activated by laser light emerged in medicine. Photosensitization is the process of initiating a reaction using a molecule called photosensitizer (PS) able to absorb a radiation and transfer this energy to reagents. Using this absorption property, numerous natural or synthetic molecules can be used as therapeutic agents from chemical reactions occurring after their activation by radiations. a. Photophysical mechanism Photodynamic activity of a PS is based on the interaction between energy delivered by the light and its absorption by the molecule (see figure 1.4). The absorption of this energy by an electron excites the PS, moving an electron to an upper level of energy. This absorption can only occur if the incident energy is exactly the same required to make this transition.

Consequently, only allowed electronic transitions can be used to activate the PS. Using the Planck's relation between energy and wavelength, these electronic transitions follow a specific absorption spectrum. Each absorption peak matches with an allowed electronic transition of each atom constituting the PS. Thus, in order to excite the PS, it is required to use the correct wavelength. There are two different processes to return to the ground state: radiative and non-radiative transitions.

On each relaxation process, dissipation of energy is required. Non-radiative transitions do not involve the light radiation emission but rather a heat emission for example.

Conversely, radiative transition emits a photon with an energy that corresponds to the gap between the two allowed energy levels. In the case of a radiative process, the emitted photon has a lower energy than the incident one due to vibrational relaxation and internal conversion.

Two different types of relaxation lead to light emission: fluorescence or phosphorescence (depending on the initial singlet or triplet state). The PS, initially at the ground and singlet state, absorbs a photon and goes straight to an excited singlet state. Using vibrational relaxation, the PS energy decreases slightly to reach the lower allowed transition of its excited state. At this stage, this singlet state can move to a triplet state. Fluorescence occurs with the singlet state relaxation, phosphorescence with the triplet one. In the case of the triplet state, photo-oxidation mechanisms produce reactive oxygen species (ROS) (such as hydrogen peroxide, hydroxyl and Superoxide radicals, singlet oxygen) and free radicals. These species lead to the degradation of several cell components and induce the cell death. Limitation appears when ROS, which is produced by the PS activation, damages the PS. This gradual degradation of the activated PS during illumination is called photobleaching [18]. This reaction induces a reduction in PS absorption and fluorescence [19][20][21][22].

Photobleaching results from PS photodegradation, which is the modification of the molecule structure (i.e., fragmentation of the sensitizer, relocalization and photomodification without rupture of the porphyrin rings) leading to the creation of non-absorbent photoproducts [23,24].

b. First and other PS generations

One of the first PS, called first generation PS, is mainly composed of hematoporphyrin and hematoporphyrin derivative (HpD). Photofrin®, the commercial name of the HpD, was one of the main PS used in several studies and clinical trials using photodynamic mechanisms [25][26][27][28].

An ideal PS should have at the same time a strong selectivity for tumor tissues and the ability to generate ROS efficiently. It should also have a stable composition and be easily produced. In addition, the PS spectrum absorption should display several high absorption coefficients in the red or near-infrared part in order to be activated as far as possible in biological tissues. Regarding the pharmacokinetics parameters, an ideal PS should have a favorable absorption, distribution, metabolism, excretion and the lowest post-treatment photosensitivity [29].

In the light of this description, this first generation of PS had several drawbacks such as a very poor selectivity and a very long duration of photosensitivity (about 30 days). In order to circumvent those issues, two generations of PS have been developed. The porphyrins group development (aminolevulinic acid, texaphyrin, or phthalocyanine) was the main interest of the second generation. Those molecules have good tissue selectivity and a modest absorption on the range 600-700nm that allows them to be activated until several millimeters in biological tissue. They own a fluorescence property (light emission around 450nm) and a correct ROS production. Some molecules are already commercialized such as the temoporfin, (commercial name Foscan®), verteporfin (commercial name Visudyne®) and the 5aminolevulinic acid (5-ALA) (commercial name Gliolan®).

Although the second PS generation shows correct results and is deeply involved in many clinical trials, the third PS generation is currently under investigation mainly in order to increase the tissue selectivity. Several technics are under development to graft other molecules on second PS generation to better target tumor molecules. Nano-particles are also highly attractive to improve the PS selectivity using encapsulation process or being used as contrast agent for imaging [30][31][32][33]. c. Case of the 5-aminolevulinic acid Because of its growing interest in oncology field using its photophysical properties, it appears essential to highlight the 5-ALA molecule. 5-ALA is a precursor drug that is converted into the actual PS, protoporphyrin IX (PpIX), by using the heme-synthesis [34] (see figure 1.5). 5-ALA administration leads to high tumor-selective PpIX accumulation in the tumor cells and its strong ability to penetrate the blood brain barrier (BBB) also contributes to the high PpIX uptake. PpIX is a photoactive substance endogenously synthesized in the metabolic pathway for heme. In a healthy cell, the maximum rate of ALA synthesis is always lower than the maximum rate of the heme synthesis. If heme is less consumed than synthesized, the ALA synthesis is naturally inhibited. Thus, neither the intermediate products As said previously, PS are molecules that use their photophysical properties to induce chemical reactions. For some decades, researchers and clinicians have tried to put this photosensitization to good use in clinical applications and particularly in oncology.

Photodynamic therapy a. History

The very first applications of a treatment using daylight (named "heliotherapy") appeared during the antiquity with Egypt, Chinese, Indian and Greek civilizations to treat skin diseases mostly [10]. However, a therapy using an exogenous substance (extern from the metabolism) reacting with photons' sun was only observed in India during the XVth century BC. Psoralen, a natural chemical compound originating from Psoralea Corylifolia stems that is a common plant in India, allowed a repigmentation of the skin (auto-immune disease leading to large white plate on body skin). At the beginning of the XXth century, Oscar Raab was the first to introduce a scientific consideration about the interaction between light and a molecule (a dye called "acridine") in order to obtain a medical treatment [10]. The term "photodynamic therapy" (PDT) appeared for the first time in 1907 in a book written by Van Tappeiner to describe the necessity of the presence of the oxygen to obtain a chemical reaction [35]. Following this publication, Van Tappeiner leaded the first clinical trial on three patients with skin cancer using PDT. An eosin solution had first been applied on treated areas and then been irradiated using sun light or an arc lamp during several weeks. He observed an improvement of patients' condition.

b. Definition

PDT is a non-thermal energy based therapy relying on light exposure after PS accumulation in the tumor cells. A PDT effect is the result of the synergy of three elements: a PS, the energy delivered through light exposure at a specific wavelength to excite the PS and the presence of oxygen. Combination of these three components leads to the formation of ROS, including singlet oxygen, which are cytotoxic molecules that damage the PS-targeted cells.

The sequence of a PDT treatment begins with the PS administration to the patient topically, orally or by intravenous injection. Once the PS goes through the blood circulation's patient, it spreads into the organism using several proteins according to its physical and chemical characteristics. It concentrates selectively using different processes (depending on whether the drug is a precursor or the PS itself) into the neoplastic tissues. In this way, a higher PS concentration (that can vary according to the tumor-selective PS uptake) is observed in the tumor cell than in the healthy one.

Under irradiation, several reactions occur between PS and the oxygen inside the cell.

In order to excite the PS, the light wavelength used must fit into one peak absorption of the PS. In most cases, a LASER light is applied (use of the monochromaticity property) as closely as possible to the therapeutic target. Light depth penetration depends on two factors: the wavelength and the size of the light source.

The higher wavelength, the deeper the light penetrates tissues (see figure 1.6a). A compromise must be found between a minimal light absorption by tissues and a maximal light absorption by the PS.

Light penetration depends also on the size of the light source (see figure 1.6b). The larger the light source, the deeper light penetrates tissues, independently from the wavelength.

However, a limit of depth penetration is reached even if the light source continues to grow. [36]).

When using a PS with a poor selectivity, light must be applied only on areas to be treated (e.g., hematoporphyrin, Photofrin® [26,27]) and organs at risk must be defined to avoid side effects. Conversely, regarding PS with a higher selectivity (e.g the temoporfin, Foscan® or protoporphyrin IX, Gliolan® [37][38][39][40]), defining safety margins to achieve a selective treatment is not needed. Indeed, if the PS uptake is low, cells will not be killed.

Several lighting modalities are available and can vary according to the topology of the target zone (see figure 1.7):

-Optical fibers, called "frontal diffusers", emit a lighting cone that can vary according the numerical aperture of the fiber. They are mostly used in topical illumination or interstitially but remain little used in clinic, -Optical fibers, called "cylindrical diffusers", allow to illuminate a larger volume than a frontal diffuser. They are used for interstitial or cavity treatments, -Lighting devices with a shape of one or several balloons (e.g., the Freiburg flap) can be used to illuminate cavity with a more homogeneous light distribution than fibers, -LED panels are often used in dermatology applications to illuminate large surfaces, -Lighting fabrics are constituted by optical fibers braided or knitted. This modality allows to obtain a more homogeneous light distribution on non-plane surfaces compared to LED panels, -The "daylight PDT" uses the light produced by the sun to bring a low power light for dermatology treatments. The creation of singlet oxygen, main cytotoxic product of a PDT reaction, induces tough modifications inside the intracellular media: lipids oxidation, spreading of amino acid and cross-linkage of protein components of cell membranes [42]. These modifications can impair permeability of cell's membrane, provoke a loss of fluidity and inhibit enzyme receptors on cell membrane, which leads to the cell death. Lifetime of the singlet oxygen is very short (approximately 50 nanoseconds). Its distance of diffusion is very short (approximately 20 nm), which restricts the PDT reaction inside cell's organelles only. Several organelles concentrate PS: mitochondria, lysosome, endoplasmic reticulum, Golgi apparatus and inside cell membrane [43,44].

Two mechanisms of cell death can occur during a PDT treatment: apoptosis and necrosis (see figure 1.8).

Apoptosis defines the programmed cell death. This type of cell death occurs naturally inside the organism and allows to evacuate damaged cells. Multiplication of damaged cells could lead to the tumor creation and development. Thus, apoptosis is one of the filters that can stop the tumor emergence and progression.

Necrosis defines the non-programmed cell death. This destruction is not determined by internal cell actions but by external cell factors leading to the cell's membrane breach.

Because necrosis produces a strong inflammatory response, apoptosis is generally preferred to decrease clinical complications.

d. Indirect mechanisms: Vascular and immune

Currently under investigation, other mechanisms with delayed actions appear also after a PDT treatment.

During a PDT treatment, a large proportion of oxygen inside tissues is consumed by the chemical reaction. The weak oxygen concentration can lead to an oxidative stress and provoke an ischemia of the treated tumor. To overcome this lack of oxygen supply and to continue its development, the tumor creates new vessels (neo-angiogenesis). Vascular endothelial growth factor (VEGF) are spread into the extracellular media during hypoxia, which can amplify the angiogenesis process [45]. An anti-angiogenesis treatment is often coupled to a PDT treatment to annihilate the process and to decrease the tumor recurrence risk. However, recent studies measure a decrease of VEGF secretion post-PDT on some cell lines [46][START_REF] Sharwani | Photodynamic therapy down-regulates the invasion promoting factors in human oral cancer[END_REF][START_REF] Kawczyk-Krupka | ALA-induced photodynamic effect on vitality, apoptosis, and secretion of vascular endothelial growth factor (VEGF) by colon cancer cells in normoxic environment in vitro[END_REF][START_REF] Xie | Effect of 5-ALA-PDT on VEGF and PCNA expression in human NPCbearing nude mice[END_REF][START_REF] Yi | Photodynamic therapy mediated by 5-aminolevulinic acid suppresses gliomas growth by decreasing the microvessels[END_REF]. These results demonstrate the significance of choosing a PS that best fits the targeted cells by PDT [42].

A growing number of studies about the long-term effect of the PDT using the activation of the immune system is observed. During the inflammatory response induced by the destruction of cells by reactive species, Damage-Associated Molecular Patterns (DAMPs) are secreted and can be recognized by Pattern Recognition Receptors (PRR). These PRR are receptors of the innate immune system. Cytokines (interleukine and chemokynes are created and differentiate immune cells to provoke an immune response to inflammation [START_REF] Maeding | Boosting Tumor-Specific Immunity Using PDT[END_REF][START_REF] Yang | Targeting Antitumor Immune Response for Enhancing the Efficacy of Photodynamic Therapy of Cancer: Recent Advances and Future Perspectives[END_REF]. Also, DAMPs activate the differentiation of T cells (CD4+ and CD8+) and B cells to initiate a response of the adaptive immune system [START_REF] Wachowska | Immunological aspects of antitumor photodynamic therapy outcome[END_REF]. Activation of the innate and adaptive immune system using inflammation due to reactive oxygen species (ROS) during PDT treatment might provoke a long-term destruction of the tumor and, consequently, decreases the recurrence risk. Figure 1.8: Illustration of cell death mechanism induced by PDT treatment (inspired by [START_REF] Pizova | Photodynamic therapy for enhancing antitumour immunity[END_REF][START_REF] Castano | Mechanisms in photodynamic therapy: Part three-Photosensitizer pharmacokinetics, biodistribution, tumor localization and modes of tumor destruction[END_REF]).

e. Applications

A growing number of scientific papers about PDT is observed: 443 results in 2000, 1425 results in 2015 (obtained with MEDLINE database (Medical Literature Analysis and Retrieval System Online) using the keywords "photodynamic therapy"). Applications are currently under investigation in the following fields:

-Oncology: head and neck cancer, pulmonology, urology, neurosurgery, gastroenterology, dermatology, gynecology [START_REF] Zhu | The role of photodynamic therapy (PDT) physics[END_REF] -Antimicrobial agent: ulcer treatment, sinusitis, acne and skin burn.

-Ophthalmology: retinoblastoma treatment

II. Glioblastoma

Gliomas are the most common primary brain tumors in adults. They are malignant tumors based on histologic similarity to mature glial cells, including astrocytes and oligodendrocytes or a mix of both cell types. According to the latest World Health Organization (WHO) classification [START_REF] Louis | The 2016 World Health Organization Classification of Tumors of the Central Nervous System: a summary[END_REF], they are categorized from Grade I to grade IV depending on their malignancy (see figure 1.9). Low-grade gliomas (LGG) include pylocytic astrocytoma (WHO grade I) and astrocytoma, oligodendroglioma (WHO Grade II). Highgrade gliomas (HGG) include anaplastic astrocytoma (WHO grade III), anaplastic oligodendroglioma (WHO grade III), and Glioblastoma (GBM, WHO Grade IV). Figure 1.9: Gliomas classification according to the world health organization.

Description and latest classification

Among them, GBM have the most dismal prognosis. GBM is a malignant astrocytic tumor with the most common glioma histology. It represents the most frequent primary brain tumor in adults [1][2][3] with an annual incidence of about 3 to 5 cases for 100,000 persons [4][5][6].

Prevalence is estimated at about 3 cases for 100,000 persons each year [START_REF] Asklund | Brain tumors in Sweden: data from a population-based registry 1999-2012[END_REF][START_REF] Brodbelt | Glioblastoma in England: 2007-2011[END_REF][START_REF] Gramatzki | Glioblastoma in the Canton of Zurich, Switzerland revisited: 2005 to 2009[END_REF][START_REF] Muquit | Socio-economic characteristics of patients with glioblastoma multiforme[END_REF][START_REF] Nava | Survival effect of first-and second-line treatments for patients with primary glioblastoma: a cohort study from a prospective registry, 1997-2010[END_REF]. Although it may occur at any age, 70% of patients are between 45 and 70 years of age [START_REF] Ostrom | CBTRUS statistical report: Primary brain and central nervous system tumors diagnosed in the United States in 2006-2010[END_REF]. GBM appears usually in the brain hemispheres, but can be found anywhere in the central nervous system.

The disease often grows rapidly (over 2 to 3 months), except when the glioblastoma develops within a pre-existing low-grade astrocytoma (secondary glioblastoma). Approximately 90% of GBM are primary tumors [START_REF] Chen | Isocitrate Dehydrogenase (IDH)1/2 Mutations as Prognostic Markers in Patients With Glioblastomas[END_REF]. GBM is characterized by its diffuse property and its strong inter and intra-tissue heterogeneity. Regarding the etiology of the disease, no factor has been reported, except the exposition to ionizing radiations in the context of brain tumor radiotherapy [2,[START_REF] Ricard | Primary brain tumours in adults[END_REF].

Prognosis is poor, especially in the absence of gross total resection, in elderly patients and in case of severe neurological impairments. Lacroix et al. published in 2001 [66] a study including 416 patients with GBM. A 53 ± 14 years mean age was reported, whom 63% of male and 37% female patients. Currently, median survival varied between 14 to 15 months according the extent of tumor resection percentage [4,[START_REF] Shah | Survival Trends in Elderly Patients with Glioblastoma in the United States: a Population-based Study[END_REF][START_REF] Johnson | Glioblastoma survival in the United States before and during the temozolomide era[END_REF][START_REF] Stupp | Radiotherapy plus concomitant and adjuvant temozolomide for glioblastoma[END_REF].

a. Tumor development

The diagnosis is most often achieved with an MRI imaging following the manifestation of patient's symptoms. A strong contrast is often observed: a necrotic core surrounded by a hyper-intense signal showing the edema and tumor infiltration [START_REF] Velazquez | Fully automatic GBM segmentation in the TCGA-GBM dataset: Prognosis and correlation with VASARI features[END_REF]. GBM distinction from the other grade III gliomas is realized using the histologic profile. Criteria include hypercellularity, nuclear atypia, and mitotic activity [START_REF] Aldape | Glioblastoma: pathology, molecular mechanisms and markers[END_REF]. Furthermore, diagnosis requires either microvascular proliferation and/or tumor necrosis. However, many aspects of these histologic features remain poorly correlated with GBM genetic markers. Main primary brain tumor markers are the loss of 10q chromosome, amplification of the gene coding for the EFGR gene and the inactivation of the CDKN2A/P16 gene. However, a mutation on the IDH (Isocitrate DeHydrogenase) gene is observed on secondary GBM. In addition, among IDH wild type high-grade, the key molecular chromosomal changes are shared between histologic grade III (anaplastic astrocytoma) and GBM (histologic grade IV) tumors. Histopathologically, several patterns exist, including giant cell GBM, small cell GBM, and gliosarcoma. Gliosarcoma cell can be observed at initial diagnosis or at recurrence, and appears to have some similar genetic aberrations as GBM. The term "GBM with oligodendroglioma component" (GBM-O) is sometimes specified where the tumor appears regionally similar to anaplastic oligodendroglioma. These tumors are easily distinguished from GBM by the presence or absence of 1p/19q co-deletion. GBM-O is distinguished from anaplastic oligodendroglioma by the absence of 1p/19q deletion, and the presence of IDH mutation and 1p/19q deletion effectively defines anaplastic oligodendroglioma and is therefore incompatible with the diagnosis of GBM-O.

b. Standard of care

The first-line treatment is usually surgery (advised by the European Society for Medical Oncology (ESMO) [START_REF] Stupp | High-grade glioma: ESMO Clinical Practice Guidelines for diagnosis, treatment and follow-up[END_REF]), either to confirm the diagnosis with a biopsy (in order to assess the grade of the glioma) or to remove as much as possible of the tumor (see figure 1.10). Resection plays a crucial role in the improvement of patient's survival. Currently, aggressive resection is encouraged to realize a maximal extent of resection (EOR) using intraoperative neuro-navigation tools. This technic enables to preserve organs at risk and functional areas to provide the patient a better quality of life. A resection surgery considered as "complete" can improve the progression free survival (PFS) of about 3 to 5 months [START_REF] Lacroix | A multivariate analysis of 416 patients with glioblastoma multiforme: prognosis, extent of resection, and survival[END_REF][START_REF] Kuhnt | Correlation of the extent of tumor volume resection and patient survival in surgery of glioblastoma multiforme with high-field intraoperative MRI guidance[END_REF][START_REF] Stummer | Extent of resection and survival in glioblastoma multiforme: identification of and adjustment for bias[END_REF].

Increase of the EOR and prolonging of the progression-free survival (PFS) have been observed when resection is carried under blue light after administration of 5-ALA (Fluorescence Guided Resection, FGR). A PFS improvement of about 2 months was observed using FGR. Because PpIX is a highly selective photosensitizer, the PS emits a fluorescence light (635 nm) if the resection cavity is exposed to blue light (375 -440 nm) [START_REF] Stepp | ALA and Malignant Glioma: Fluorescence-Guided Resection and Photodynamic Treatment[END_REF][START_REF] Leroy | Fluorescence guided resection and glioblastoma in 2015: A review[END_REF]. This light emission helps the surgeon detecting remaining tumor tissue at surface of the resection cavity and therefore, to improve the EOR. Unfortunately, improvement of overall survival cannot be currently statistically proved although there is a trend toward increase of survival [START_REF] Guyotat | 5-Aminolevulinic Acid-Protoporphyrin IX Fluorescence-Guided Surgery of High-Grade Gliomas: A Systematic Review[END_REF].

Once the surgery or biopsy realized, adjuvant treatments are administered according to the grade and the genetic profile of the tumor with different treatment options (figure 1.10): The loss of heterozygosity of the 1p19q gene [START_REF] Gadji | Is the 1p/19q deletion a diagnostic marker of oligodendrogliomas?[END_REF][START_REF] Zhao | Loss of heterozygosity 1p/19q and survival in glioma: a meta-analysis[END_REF], IDH mutation and methylguanine methyltransferase (MGMT) methylation enzyme (enzyme used in DNA reparation process) are the three major keys points occurring into the decision algorithm of adjuvant treatments.

In terms of survival, the benefits from adjuvant treatments after resection surgery are significant, although they remain modest. Because tumor cells infiltrate the whole brain, a full resection cannot be performed. Radiotherapy combined with chemotherapy help to complete the treatment.

In case of relapse, the WHO does not establish a standard of care. Generally, a secondline of chemotherapy or surgery is achieved. c. Limits of the standard treatment and alternative treatments Although standard treatment for patients with GBM, that includes surgical resection, radio and chemotherapy, allows to improve the patient's survival, patients cannot be cured of their disease. These modest results are in part due to the cellular composition and more particularly, to the presence of self-renewing, tumorigenic stem cells that contribute to tumor initiation and recurrence [START_REF] Lathia | Cancer stem cells in glioblastoma[END_REF][START_REF] Swetha | The role of cancer stem cells in glioblastoma[END_REF]. Furthermore, conventional chemotherapies fail to bring an effective response to GBM treatment due to a poor BBB penetration. In this context, several alternative therapies appear using different treatment ways [START_REF] Thomas | Emerging therapies for glioblastoma[END_REF]:

-Immunotherapy: this recent treatment modality is based on the immune system activation. In GBM treatment context, checkpoint molecules and monoclonal antibodies are targeted. Checkpoint molecules inhibition is aimed at enhancing immune activation. Monoclonal antibodies target tumor-specific antigens to cause tumor cell death [START_REF] Cohen-Inbar | Focused ultrasound-aided immunomodulation in glioblastoma multiforme: a therapeutic concept[END_REF][START_REF] Hung | Biomarkers and Immunotherapeutic Targets in Glioblastoma[END_REF].

-Epigenetic: this therapy affects the genes expression. In GBM treatment context, Histone deacetylase (HDAC) inhibitors appear as a promising way to stop tumorigenesis [START_REF] Lee | Advances in epigenetic glioblastoma therapy[END_REF]. HDAC expression is related directly to the patient's survival.

HDAC inhibition can provoke both direct effects (apoptosis, enhanced ROS production) and indirect effects (angiogenesis, metastasis and glucose utilization inhibition).

-HIFU: the therapy named High-Intensity Focused Ultrasound is mostly combined with chemotherapy for the GBM treatment. This therapy is based on the emission of multiple ultrasound (approximately 1 Mhz) beams focused on the same spot in order to create a thermal and/or a mechanical effect [START_REF] Cohen-Inbar | Focused ultrasound-aided immunomodulation in glioblastoma multiforme: a therapeutic concept[END_REF][START_REF] Yang | Chemotherapy of glioblastoma by targeted liposomal platinum compounds with focused ultrasound[END_REF][START_REF] Phenix | High intensity focused ultrasound technology, its scope and applications in therapy and drug delivery[END_REF][START_REF] Yang | Treating glioblastoma multiforme with selective high-dose liposomal doxorubicin chemotherapy induced by repeated focused ultrasound[END_REF][START_REF] Yang | Pulsed High-Intensity Focused Ultrasound Enhances the Relative Permeability of the Blood-Tumor Barrier in a Glioma-Bearing Rat Model[END_REF]. The main goal of this approach is to disrupt the BBB to allow drug used in chemotherapies to penetrate into tumor tissues.

-Proton therapy: currently, conventional radiotherapies do not deposit a precise enough dose, as brain tumors treatment requires. However, proton therapy brings a successful way to obtain a high dose in a very selective volume and to reduce radio-toxicities [START_REF] Von Neubeck | Glioblastoma multiforme: emerging treatments and stratification markers beyond new drugs[END_REF]. However, like photon radiotherapy, proton therapy can target only the visible part of the tumor with a small margin of macroscopic healthy brain tumor.

III. PDT treatment of HGG 1. History

To our knowledge, PDT to treat brain tumors was first studied in 1972 [11]. From the observation of Haematoporphyrin concentration in neoplasm tissues, the authors grafted 20 rats with GBM cells. Tumor-selective PS uptake along with the death of tumor cells after light exposure was observed. However, the authors failed to obtain total tumor destruction. They emitted the hypothesis that the inhomogeneous PS or the light distribution could lead to this incomplete treatment.

In 1980, the first clinical trial using intraoperative PDT to treat HGG was realized by Perria et al. [12]. Nine patients were treated using a He-Ne laser (632.8 nm) to illuminate the resection cavity with hematoporphyrin derivative (HpD) PS. Each patient received a fluence of 9 J/cm² with an injection between 2.5 to 10.0 mg/kg of HpD. No adjuvant treatments (radio or chemotherapy) were delivered to the patients. However, survival remains very weak: between 6 to 44 weeks.

One year later, in 1981, a report on the first phase I clinical trial using interstitial PDT treatment was published by Laws et al.. Five patients were treated with 5 mg/kg of HpD and were protected from direct or indirect sunlight during at least 2 weeks [START_REF] Laws | Photoradiation Therapy in the Treatment of Malignant Brain Tumors: A Phase I (Feasibility) Study[END_REF]. A single optical fiber (400µm diameter with a numerical aperture of 0.35) was inserted at the center of the tumor using stereotactic localization. The argon-pumped laser provided red light (630 nm) with a power of 300 to 400mW delivered at the end of the fiber. Authors administered a 300mW red light during 45 minutes. Three grade III and two grade IV were treated. The authors concluded that a CT-guided needles insertion might be interesting during treatment planning to insert multiple fibers in order to reach an optimal tumor exposure to laser light.

Reaction to PDT could be observed in CT post-treatment.

Interstitial PDT

When the surgical resection cannot be achieved because of surgical constraints (e.g., functional area, relapse), PDT delivered interstitially is investigated (approximately 15% of all cases [START_REF] Fazeny-Dorner | Survival and prognostic factors of patients with unresectable glioblastoma multiforme[END_REF][START_REF] Paquette | Prognostic value of health-related quality of life for death risk stratification in patients with unresectable glioblastoma[END_REF][START_REF] Mirimanoff | Radiotherapy and temozolomide for newly diagnosed glioblastoma: recursive partitioning analysis of the EORTC 26981/22981-NCIC CE3 phase III randomized trial[END_REF]). This modality plans to insert optical fibers under stereotactic conditions and to vehicle light directly into the tumor (see figure 1.11a).

a. Instrumentation

In most studies, optical fibers are used to deliver laser light through brain tissues. Two types of medical devices are used: frontal and cylindrical diffusers as mentioned before (see figure 1.7). They are mainly made with silica core (0.4 to 0.6mm diameter) surrounded by a plastic casing. Frontal diffusers (or bare-end fibers) are characterized by their numerical aperture, which can vary from 0.2 to 0.4. This feature defines the size of the cone beam: the higher the numerical aperture, the higher the angle of the cone is. This type of diffuser is mostly used in dermatology [START_REF] Johansson | In vivo measurement of parameters of dosimetric importance during interstitial photodynamic therapy of thick skin tumors[END_REF], urology [START_REF] Johansson | Realtime light dosimetry software tools for interstitial photodynamic therapy of the human prostate[END_REF], gastroenterology [START_REF] Wolfsen | Bare fiber photodynamic therapy using porfimer sodium for esophageal disease[END_REF], pre-clinical trials on small animals [START_REF] Tetard | Interstitial 5-ALA photodynamic therapy and glioblastoma: preclinical model development and preliminary results[END_REF] or in optical measurements [START_REF] Beck | Clinical Determination of Tissue Optical Properties in vivo by Spatially Resolved Reflectance Measurements[END_REF]. Cylindrical diffusers produce an ellipsoid illumination around the extremity composed of a diffusing tip. Currently, this type of diffuser is more used in clinical applications because it maximizes the irradiated volume during the procedure and thus, minimizes both the number of fiber to insert and the treatment time [100].

With technological improvements (imaging and computing performance), a few treatment planning systems (TPS) have been developed in PDT field: pleural application Improvement of accuracy of surgical gesture in neurosurgery is mainly due to the apparition of the stereotactic surgery method [107]. In 1908, Horsley and Clarke published a paper describing the first apparatus dedicated to a monkey head and allowing to localize a point in a 3D Cartesian coordinate system. The first trial of stereotactic devices on humans was realized by Speigel et al. in 1947 [108]. In 1951, Lars Leksell described an apparatus setup on the patient's head, defining a 3D polar coordinate system much easier to use during surgical procedure [109]. At the same period, Jean Talairach and Pierre Tournoux developed also a stereotaxic frame [110,111] to create a brain atlas [112]. Since their respective commercialization, several other systems have been proposed to improve the technic and adapt it to the advent of imaging system (CT and MRI). This tool could be a real advantage in interstitial PDT treatment for both planning step and surgical gesture and, furthermore, for computer-assisted which is on the rise. Robot-assisted stereotactic neurosurgery is an technology allowing high accuracy treatment on laser interstitial thermal therapy, responsive neuro-stimulation or stereo electroencephalography [113,114]. Replacing conventional tools (e.g., electrodes, needles) by optical fibers, an accurate interstitial PDT treatment could be delivered by such robot-aided system.

b. Clinical trials

The following table represents the major clinical trials of interstitial PDT treatment of GBM.

Authors

Intraoperative PDT

Intraoperative PDT treatment has been the first PDT modality used to treat brain tumors (see figure 1.11b). This lighting modality aims at illuminating the cavity created by the surgical resection in order to treat the remaining tumor tissue that cannot be removed. It also aims at treating the healthy brain invaded by infiltrating tumor cells inevitably.

a. Instrumentation

The first clinical trial of brain tumors using PDT achieved by Perria et al. in 1980 had been performed during the surgical procedure [12]. A He-Ne Laser (632. Some clinical trials used a single optical fiber placed at the center of the cavity.

Differences appear between frontal diffuser, cylindrical diffuser and spherical diffuser. In many cases, the resection cavity is filled with diffusing solution.

Another type of irradiator was proposed by Schmidt et al. in 1996 using properties of LED [123]. An inflatable balloon-shape device was designed to fit into the resection cavity to deliver light provided by a light source composed of 144 LED chips. These chips were arranged as a cylinder and placed at the center of the balloon. Three inputs were created in the tube fixed to the balloon: firstly to provide electricity for the LED tip, secondly to fill the balloon with diffusing liquid and thirdly, to cool the LED tips using sterile water. A total light power of 1.0W was measured with a peak emission at 677nm and with a wide bandwidth of about 20 -25nm. This device was evaluated in clinical trial and compared to laser beam.

Recently

Chapter conclusion

With the improvement and miniaturization of LASER device, the use of this technology increases in the medical field whether it is in esthetic surgery, ophthalmology, dermatology than cancerology.

Despite strong efforts accomplished to improve GBM care, most of studies revealed that the current standard of care does not bring an effective solution for the recovery of the patient. Moreover, adjuvant treatments such as radio and chemotherapy have a modest impact on the extension of the PFS and even more on the overall survival.

Among recent therapies, the use of LASER with the PDT has shown an interest and significance for the management of GBM, even when the surgical resection cannot be achieved. Several studies have already been publishing with encouraging results. However, PDT has not been adapted to clinical routine for several issues yet. No consensus has been defined on the use of a photosensitizer dedicated to GBM treatment. A lack of reliable and reproducible therapy schemes regarding the clinical devices, light delivery systems and TPSs is also observed between clinical studies. All of these disparities cause an absence of randomized and controlled multicentre clinical trials on the GBM 5-ALA PDT treatment.

Thus, the main objective of this thesis work is to provide common tools (software and light delivery systems) for PDT in the GBM management.

Chapter II Photodynamic therapy dosimetry Predictive dosimetry models the three major features PDT components (light propagation, photosensitizer distribution and oxygen concentration) in order to plan the PDT treatment delivery. The monitoring of PDT treatment aims to evaluate the preoperative planning and if necessary, to optimize intraoperatively the PDT treatment.

To obtain a relevant dosimetry, a comprehensive understanding of light interaction within biological media is essential. Thus, physics interactions between light and matter are presented from microscopic to macroscopic scale. Since several models have been developed to compute these dosimetric parameters, a state of the art with a focus on light propagation computing through biological media is introduced. Additionally, with improvement of measuring systems, monitoring markers used in PDT monitoring are shortly described.

I. Light and matter

In this first part, only the particle theory of the wave-particle duality is considered to describe processes occurring during photons history. Throughout this manuscript, "light" and "photon" refer to theories applied in the range of non-ionizing radiations energies.

Optical properties

From a simplified point of view, light constituted of photons may be subjected to two interactions during its history: it can be absorbed or not. When not absorbed, the photon may be deviated from its initial direction and takes another direction, and even turns back toward its birthplace.

Several parameters describe and quantify the interactions of photons within matter.

They define matter properties and how they perturb the photon propagation. It appears necessary to define them for the future developments in this thesis.

a. Absorption

The first parameter is the absorption coefficient (µa) [m -1 ] of a given media. Let us consider the molecules, called chromophores, composing this media. A surface called effective cross-section σa [m 2 ], contained into these chromophores, absorbs light (see figure 2.1). The volume density of the medium containing the chromophores is ρa [m -3 ] (i.e number of chromophores per volume of media). The absorption coefficient is thus defined by:

µ 𝑎 = ρ 𝑎 . σ 𝑎 (2.1)
Furthermore, the probability Psurvival of photon's survival after crossing a length L [m] within the media is:

𝑃 𝑠𝑢𝑟𝑣𝑖𝑣𝑎𝑙 = 𝑒 -µ 𝑎 .𝐿 (2.2)

b. Scattering

Light may be absorbed or deviated from its initial direction. The coefficient describing this phenomenon is named scattering coefficient. This phenomenon results from the interaction of a photon with an obstacle on its pathway. Mostly, the photon wavelength is not affected during this process. The term "elastic" is employed to describe this collision during which the total kinetic energy of the system photoncollided particle is conserved. This system can be considered as a harmonic oscillator model, i.e., two masses tied by a spring (see figure 2.2). When the incident photon wavelength matches exactly with the resonant frequency of the harmonic oscillator system fresonance, the entire energy of the incident photon is absorbed by the system. Also, if the incident photon wavelength does not match exactly with the resonant frequency of the harmonic oscillator system, the system enters in a forced oscillation leading to a very poor energy transfer. Nevertheless, an oscillating system reacts as an antenna and radiates its energy outward. Consequently, the incident photon can provoke oscillation of the mass-spring system, which remits energy in a different direction from the incident direction; the photon is scattered.

If the size of the obstacle is close to the photon wavelength, this one has a higher probability to be scattered. Two theories can be applied according to the size of the encountered structure. If the size of this latter is of the same scale as the incident photon wavelength, the Mie Theory describes the scattering process. If the size of the encountered structure is smaller than the length of the incident photon wavelength, the Rayleigh scattering theory describes scattering process (see figure 2.3). The Mie theory is a particular solution of Maxwell's equations used to describe the elastic scattering. As for Rayleigh scattering, it approximates the limit conditions of the Mie theory, especially for particles much smaller than the incident radiation wavelength. Regarding the wavelength scale mostly used in PDT (between 400 to 900 nm), the Mie theory is applied to explain the photon's behavior during the scattering.

The parameter that describes the scattering event is the scattering coefficient (µs) [m

µ 𝑠 = ρ 𝑠 . σ 𝑠 (2.3)
Furthermore, the probability Ptransmission of photon crossing a length L [m] of this media without scattering event is: The total attenuation coefficient µt is the sum of the absorption and scattering coefficients:

𝑃 𝑡𝑟𝑎𝑛𝑠𝑚𝑖𝑠𝑠𝑖𝑜𝑛 = 𝑒 -µ 𝑠 .𝐿 (2.4)
µ 𝑡 = µ 𝑎 + µ 𝑠 (2.5)
Finally, on an infinitesimal path length ΔL, the probability of a photon to interact (absorption or scattering process) is given by:

𝑃 𝑖𝑛𝑡𝑒𝑟𝑎𝑐𝑡 = 𝑒 -µ 𝑡 .ΔL (2.6) c. Anisotropy
After a scattering event, the photon has a new trajectory. The coefficient of anisotropy g [dimensionless] describes the behavior of the photon's direction, and more specifically the angle between the incident trajectory and the new one. Let us consider a photon scattered by a structure. The angle θ corresponds to the difference between the initial trajectory and the deviated one. The anisotropic coefficient defines the mean cosines value of the expected angle θ after a scattering event:

𝑔 = < 𝑐𝑜𝑠 (𝜃) > (2.7)
An ideal isotropic scattering process is described by g = 0 and a total forward scattering process is described by g = 1. The scattering phase function p(θ) [sr -1 ] designates the photon probability to be scattered into a solid angle oriented at the angle θ according to the value of the anisotropy g (see figure 2 The anisotropic coefficient can be described as [135]: 2.9) This equation better represents the notion of probability (due to the scattering phase function)

𝑔 = ∫ 𝑝( 𝜋 0 𝜃). cos(𝜃). 2𝜋. sin(𝜃) . 𝑑𝜃 (
in the anisotropic coefficient (see figure 2.5).

The reduced scattering coefficient µ's [m -1 ] is an important factor often used in optic to take into account the anisotropy of the media:

µ′ 𝑠 = µ 𝑠 . (1 -𝑔) (2.10)
The inverse of the reduced scattering coefficient describes the mean distance of a photon according to the anisotropy of the medium. This shortcut, called reduced mean free path, is equal to several infinitesimal photon displacements in a high scattering media (the absorption coefficient must be negligible when compared to the scattering coefficient) (see figure 2.6).

Radiometric parameters

In this part, non-exhaustive definitions of radiometric parameters are given. Only parameters mostly used in PDT dosimetry were selected to be highlighted.

a. Light sources

Light sources are mostly defined according to their illumination type (e.g., in watt for frontal diffusors, in watt per centimeter for cylindrical diffusors). This trivial definition of the light source can be detailed in several other concepts used in biomedical optics.

The photon density N [photons.m -3 ] emitted from this light source can be established.

It characterizes the number of photons inside a unit volume (defined by a solid angle dω directed by the unit vector û).

Using this photon density, the radiance term L [W.m -2 .sr -1 ] introduces the notion of energy in this previous density. It represents the power per unit area (and not per unit volume)

of photons crossing this area, which is oriented perpendicularly to the unit vector 𝑠̂ and delimited by the solid angle dω (see figure 2.8). "For the novice and even for many in the field, it goes against intuition that the fluence rate, φ [W/m 2 ] within tissue can be larger than the irradiance E0 [W/m 2 ]."

The irradiance defines the power incident delivered by a light beam on an area situated at the targeted tissue surface only, similarly to a light beam hitting a screen. This expression considers only the geometric configuration of the light source and the illuminated surface and is mainly employed in dermatology because most of illumination modalities are delivered topically. The amount of energy delivered per second at the surface of the skin can be expressed as an irradiance term. Nevertheless, in many biomedical optics applications, the fluence rate has to be employed.

Fluence rate considers the radiance obtained over all directions. It expresses the power incident delivered by a light beam inside a small sphere, on the cross-sectional area of this small sphere. Thus, the area can be situated not only on the surface but also inside tissues.

Due to scattering events occurring in highly diffusing medium, one photon can cross several times the cross-sectional area where the fluence rate is estimated. This difference explains values of the fluence rate that can be higher than the irradiance, even if they are expressed both in the same unit [W.m -2 ] [136, 137] (see figure 2.9). The specular reflectance can be estimated using the Fresnel equation:

𝑅 𝑠 = 1 2 . [ 𝑡𝑎𝑛²(𝜃 1 -𝜃 2 ) 𝑡𝑎𝑛²(𝜃 1 + 𝜃 2 ) + 𝑠𝑖𝑛²(𝜃 1 -𝜃 2 ) 𝑠𝑖𝑛²(𝜃 1 + 𝜃 2 ) ] with 𝜃 1 ≠ 0 (2.16)
When the incident angle θ1 is null, i.e., the incident angle is perpendicular to the interface, the expression becomes:

𝑅 𝑠 = (𝑛 1 -𝑛 2 )² (𝑛 1 + 𝑛 2 )² with 𝜃 1 = 0 (2.17) Light transmission T [dimensionless]
is the ratio of light that is not reflected at the interface and penetrates the second media, can be estimated using the specular reflectance:

𝑇 = 1 -𝑅 𝑠 (2.18)
Derived from the light transmission, the transmittance is the proportion of light absorbed by a medium sample with a known thickness.

II. PDT dosimetry: state of the art

In their review, Pogue et al. defines the dosimetry as [138]:

"The quantitative planning and verification processes used to prescribe and verify a patient treatment, which will ensure that the targeted areas are given appropriate dose to kill the tissue, and to prevent over-treatment to normal tissues, avoiding morbidity issues."

In a PDT context, this dosimetry definition describes both the planning of the treatment delivery (or predictive dosimetry) and the monitoring during the treatment.

Planning includes prediction of light propagation, PS concentration and tissue oxygen concentration prior to the treatment delivery (relative to patient treatment prescription).

Monitoring is thus the verification of preoperative planning (or prescription). In this second paragraph, methods applied to perform the PDT dosimetry, both for planning and monitoring are listed. Regarding the planning, modeling of the light delivery, photosensitizer distribution and oxygen concentration are presented with a focus on light propagation computing. Then, monitoring approaches such as fluorescence and singlet oxygen detection are detailed in the last section.

Predictive dosimetry

Predictive dosimetry is composed of three features that are relative to the three major PDT components: the light propagation, the photosensitizer distribution and the oxygen concentration. Modeling these three parameters enables to plan the PDT treatment delivery.

a. Light delivery computing

Three different models may be used to compute the light distribution in steady-state inside biological media. The analytical models enable to estimate quickly the fluence rate distribution in a homogeneous medium. Numerical models are most often dedicated to compute the light propagation in more complex tissues. Monte-Carlo is a physical model used to compute light transport in heterogeneous media and can be consider as the reference model.

i. Analytical models

The analytical models allow estimating the fluence rate produced by a light source using mathematical equations. All developed models originate from the Boltzmann transport equation [139]:

1 𝑐 . 𝜕𝐿(𝑟, 𝑠, 𝑡) 𝜕𝑡 + 𝛻. 𝐿(𝑟, 𝑠, 𝑡)𝑠= -(µ 𝑎 + µ 𝑠 )𝐿(𝑟, 𝑠, 𝑡) + µ 𝑠 ∬ 𝐿(𝑟, 𝑠′, 𝑡)𝑓(𝑠̂• 𝑠̂′)𝑑𝛺 ′ 4𝜋 + 𝑄(𝑟, 𝑠′, 𝑡) (2.19)
where c is the light velocity

[m.s -1 ], t is the time [s], 𝐿(𝑟, 𝑠, 𝑡) is the radiance [W.m -2 .sr -1 ],
𝑄(𝑟, 𝑠′, 𝑡) defines the power centered on the unit vector 𝑠̂ at the r position and the normalized differential scattering cross-section ∬ 𝑓(𝑠̂• 𝑠̂′)𝑑𝛺 ′ 4𝜋 = 1

The Boltzmann transport equation (equation 2.19) can be simplified using the definition of the source 𝑆(𝑟, 𝑡), the fluence rate 𝜙(𝑟, 𝑡) and the flux 𝑗(𝑟, 𝑡):

𝑆(𝑟, 𝑡) ≡ ∬ 𝑄(𝑟, 𝑠, 𝑡)𝑑𝛺 4𝜋 𝜙(𝑟, 𝑡) ≡ ∬ 𝐿(𝑟, 𝑠, 𝑡)𝑑𝛺 4𝜋 𝑗(𝑟, 𝑡) ≡ ∬ 𝐿(𝑟, 𝑠, 𝑡)𝑠d𝛺 4𝜋 (2.20)
Using the Fick's law of diffusion, the flux can be written as:

𝑗(𝑟, 𝑡) = -𝐷𝛻𝜙(𝑟, 𝑡) Where 𝐷 = 1 3[(1-𝑔)µ 𝑠 +µ 𝑎 ] = 1 3[µ′ 𝑠 +µ 𝑎 ] = 1 3µ 𝑡𝑟 (2.21)
Assuming the source as isotropic and in a steady state, the combination of these different expressions yields to the general diffusion equation in turbid medium with scattering property stronger than absorption:

1 𝑐 . 𝜕𝜙(𝑟, 𝑡) 𝜕𝑡 = 𝐷𝛻 2 𝜙(𝑟, 𝑡) -µ 𝑎 (𝑟)𝜙(𝑟, 𝑡) + 𝑆 0 (𝑟, 𝑡) (2.22)
where c is the light velocity [m.s ). This equation describes the wave propagation in space and time, emitted from a source 𝑆 0 (𝑟, 𝑡).

-1 ], 𝜙(𝑟, 𝑡) is the fluence rate [W.m -2 ], t is the time [s], D is the diffusion coefficient [m], µ 𝑎 (𝑟) is the absorption coefficient [m -1 ]
In a PDT dosimetry context, the fluence rate is most often the main variable studied in order to obtain the fluence [J.m -2 ] during the planed treatment. Consequently, only stationary problems are considered and are independent from the time t. The general diffusion equation can be simplified by dropping the t parameter and considering 𝜕𝜙(𝑟,𝑡) 𝜕𝑡 equal to zero. This leads to the following equation:

𝑆 0 (𝑟) = µ 𝑎 (𝑟)𝜙(𝑟) -𝐷𝛻 2 𝜙(𝑟) (2.23)
which is the partial differential Helmholtz's equation [140,141].

In stationary state and with a sufficient distance from light sources (considering superior to the reduced mean free path) or interfaces between several media, the Helmholtz equation can be simplified in:

𝐷𝛻 2 𝜙(𝑟) = µ 𝑎 (𝑟)𝜙(𝑟) (2.24)
Using this equation (equation 2.24), several simple light sources geometries can be modeled.

The easier model is the point source of power P with r the distance between the source and the position of the point where the fluence rate 𝜙(𝑟) is estimated:

𝜙(𝑟) = 𝜙(0) 4𝜋. 𝐷. 𝑟 . 𝑒 -𝑟 √ 𝐷 µ 𝑎 ⁄ ≈ 𝑃. 3(µ 𝑎 + µ′ 𝑠 ) 4𝜋. 𝑟 . 𝑒 -𝑟.µ 𝑒𝑓𝑓 (2.25)
More complex geometries can be modeled such as cylindrical light sources, i.e., cylindrical diffusor used in most interstitial PDT treatment. However, two equations can be found in the literature to estimate the fluence rate emitted from a light source of a length l in turbid media (equations 2.26 and 2.27). The first equation (equation 2.26) discretizes the diffusing part of the optical fiber as a sum of several point light sources [142][143][144]. Thus, calculation of the fluence rate at a distance r from the fiber is the sum of each light source contribution:

𝜙(𝑟) = ∑ P. 3. µ s ′ 4. π. 𝑟 n . e -r n .µ eff n sources n = 1 (2.26)
The second equation (equation 2.27) considers the whole fiber as a finite line light source with 2D cylindrical light emission characteristics [145,146]. It computes the fluence rate values using the minimal distance r from the fiber: (2.28)

ϕ(𝑟) = 𝑃. √ 2.
These equations allow to estimate light distribution of various light sources geometries during a very short computing time (see figure 2.11). It represents definitely the simplest approach to evaluate the light propagation and can be easily integrated in a TPS dedicated to PDT.

Nonetheless, it remains difficult to take into account heterogeneities of tissues. Fluence rate estimation closed to light source (under the reduced mean free path) and closed to boundaries requires complex mathematical approximations that make the usage of analytical models not relevant in these cases. ii. Finite elements Numerical analysis has been developed in order to resolve mathematical equations such as partial differential equation (PDE). This approach called "finite element" is largely used in engineering and physics, and subsequently in biomedical optics. Finite element can be used to solve PDE, such as the Helmholtz equation (equation 2.23). This method discretizes the space into a mesh in order to set boundaries conditions and enables to obtain numerical approximations. Currently, meshes are adapted to an expected resolution: a fine resolution is generally applied near to complex shape or close to light sources. Conversely, a coarser resolution is defined at a distance from complex geometries in order to reduce computing time. From this mesh, a resolving of diffusion equation is performed at each node constituting the mesh. Since it enables an easier modeling of shapes and it is implemented in several major modeling software such as Comsol MultiPhysics® (COMSOL Inc., Burlington, USA), only the second approach called finite element method is described in this section. A full mathematical explanation of other methods would be out of the scope of this work. Several publications already describe in more details these methods [147][148][149][150][151][152].

The From this simple consideration, a more complex model can be elaborated:

〈-𝛻 • 𝐷𝛻𝛷 𝜑 𝑖 〉 + 〈µ 𝑎 . 𝛷 𝜑 𝑖 〉 = 〈𝑆 𝜑 𝑖 〉 (2.30)
where the notation 〈𝐻〉 designates the integration of the entire field in Hilbert spaces (space with finite or infinite dimensions) and 𝛷 𝜑 𝑖 is the discrete fluence evaluated at each node i weights 𝜑 𝑖 . It does not appear necessary here to describe more in details the numerous mathematical steps to obtain the following equation (already described in [160][161][162][163][164]):

𝛷 = [𝐴] -1 𝑆 (2.31)
where A is the matrix of the diffusion connectivity between basis functions and S is the vector of sources at each node. Interpolation between nodes is performed in order to produce a continuous matrix of the diffusion equation solving. iii

. Monte-Carlo

The third model used to estimate the fluence rate in biological media is called Monte-Carlo method. This algorithm model is one of the most popular estimation algorithm and widely used in physics applications and particularly in particle physics simulations [165].

Thanks to its accuracy, this method is now in use in most TPS for planning high-energy photon radiation therapies. A well-known example to illustrate the Monte-Carlo method and understand the use of random numbers is the determination to the area of a lake 𝐴 𝑙𝑎𝑘𝑒 situated in square site 𝐴 𝑠𝑖𝑡𝑒 of known dimensions. A cannon fires 𝑥 𝑟𝑎𝑛𝑑𝑜𝑚 shots randomly in this square area. The number of balls remaining on the ground 𝑛 𝑔𝑟𝑜𝑢𝑛𝑑 is determined. Thus, the number of balls fallen in the lake are 𝑥 𝑟𝑎𝑛𝑑𝑜𝑚 -𝑛 𝑔𝑟𝑜𝑢𝑛𝑑 . The area of the lake 𝐴 𝑙𝑎𝑘𝑒 can be estimated from these numbers:

𝐴 𝑙𝑎𝑘𝑒 = 𝑥 𝑟𝑎𝑛𝑑𝑜𝑚 -𝑛 𝑔𝑟𝑜𝑢𝑛𝑑 𝑥 𝑟𝑎𝑛𝑑𝑜𝑚 . 𝐴 𝑠𝑖𝑡𝑒 (2.32)
The estimation accuracy of the area of the lake 𝐴 𝑙𝑎𝑘𝑒 improves by increasing the number of cannon fires 𝑥 𝑟𝑎𝑛𝑑𝑜𝑚 (see figure 2.15). Similar approach can be used to model light propagation in matter. Indeed, using probabilities of interaction of photons with matter previously detailed, for a given photon, each event on its pathway can be randomly defined according to these probabilities. To model the propagation of light in matter, each photon pathway is computed from birth to death.

Photons' life (creation, transport and absorption) is called history and each step of their history is recorded. Random numbers are injected many times in the program to evaluate all probability distribution functions associated to all steps of photons' history. Accordingly, by transposing this random technic to the light propagation through matter, a large number of photons had to be randomly generated in order to tend toward a good estimation of the light absorption. Using all data collection of several millions of histories, light propagation through matter can be estimated.

The first part of the third chapter is dedicated to the description of the Monte-Carlo model, the standard algorithm and all computing improvements.

b. Photosensitizer distribution

The PS uptake in tissues remains the most crucial point for efficient PDT treatment.

The modeling of the PS distribution into tissues before the PDT treatment is mainly governed by mathematical equations requiring precise knowledge of the PS pharmacokinetics [170,171]. The rate of diffusion of a PS is determined by the concentration gradients and the diffusivity. This constant depends on the tissues properties and also on the chemical properties of the PS used in the treatment [172]. The Fick's first law describes the diffusion of matter in steady state:

𝐽 = -𝐷. 𝑑𝑀 𝑑𝑥 (2.33)
Where J is the flux vector indicating the direction of the PS, D is the diffusion coefficient through mater, M is the precursor drug (e.g., 5-ALA for PpIX) concentration and x is the depth in the media.

From this simple equation, more complex equations can derivate to describe the temporal evolution of the PS and its precursor [173] and lead to the concentration of photoactive compound accumulated in the target tissue.

Although a profound comprehension of the PS uptake within tissues remains essential, according to selectivity and specificity of a given PS, the modeling of its distribution in a dosimetry context may be secondary.

c. Singlet oxygen concentration

Initially, oxygen is a prerequisite to obtain a PDT response during the illumination of targeted areas. However, oxygen level may vary spatially and dynamically due to its consumption during a PDT treatment [174][175][176][177].

Singlet oxygen is the main cytotoxic agent provoking considerable damages in tumor cells and represents a particular lever in PDT treatment efficacy. Thus, a strong effort has been made to develop a model able to predict the singlet oxygen concentration in tissues [178].

Because singlet oxygen distribution modeling depends on the reaction type implemented, most of developments have been concentrated on the type II reaction. Singlet oxygen creation requires several previous photochemical reactions and all species can be described by differential equations [174, 179] [180, 181]:

𝑑[𝑆 0 ] 𝑑𝑡 = -𝑘 0 [𝑆 0 ] -𝑘 1 [ 𝑂 1 2 ]([𝑆 0 ] + 𝛿) + 𝑘 2 [𝑇][ 𝑂 3 2 ] + 𝑘 3 [𝑆 1 ] + 𝑘 4 [𝑇] 𝑑[𝑆 1 ] 𝑑𝑡 = -(𝑘 3 + 𝑘 5 )[𝑆 1 ] + 𝑘 0 [𝑆 0 ] 𝑑[𝑇] 𝑑𝑡 = -𝑘 2 [𝑇][ 𝑂 3 2 ] -𝑘 4 [𝑇] + 𝑘 5 [𝑆 1 ] 𝑑[ 𝑂 3 2 ] 𝑑𝑡 = -𝑆 Δ 𝑘 2 [𝑇][ 𝑂 3 2 ] + 𝑘 6 [ 𝑂 1 2 ] + 𝑃 𝑑[ 𝑂 1 2 ] 𝑑𝑡 = -𝑘 1 ([𝑆 0 ] + 𝛿)[ 𝑂 1 2 ] + 𝑆 Δ 𝑘 2 [𝑇][ 𝑂 3 2 ] -𝑘 6 [ 𝑂 1 2 ] -𝑘 7 [𝐴][ 𝑂 1 2 ] 𝑑[𝐴] 𝑑𝑡 = -𝑘 7 [𝐴][ 𝑂 1 2 ]
(2.34)

where

Monitoring of the PDT treatment

Monitoring enables the evaluation of the treatment response during the illumination and the adaptation of the light delivery. With a strong uptake of tumor cells, fluorescence of the PS is the most intuitive marker to monitor the treatment response. Also, although it is still difficult to record, the monitoring of singlet oxygen luminescence is the most direct marker of the PDT dose.

a. Fluorescence detection

The emergence and now, the common practice of the fluorescence-guided resection (FGR) in clinical routine is the first application of PS fluorescence property [START_REF] Leroy | Fluorescence guided resection and glioblastoma in 2015: A review[END_REF][185][186][187][188][189].

The use of 5-ALA (commercial name Gliolan®) as contrast agent in surgery of malignant glioma has been approved for human use in Europe, Asia, Australia and USA 1 .

Thus, the measurement of the PS concentration before, during and after a PDT treatment can be included as an indicator in a clinical trial. This feature provides from direct in situ measurements of the fluorescence emitted by the PS [190][191][192][193]. Quantifying PS fluorescence has the potential to optimize PDT treatment parameters and also to monitor the PS concentration during treatment. According to PS used in the PDT treatment, high variability of concentration can appear between patients and tissues within the same tumor. Currently, spectroscopy technic is used to detect the PS fluorescence of the treated area. Several research teams have already developed several technological tools to detect this fluorescence such as fiber probe systems for deep PS distribution measurements [191,[194][195][196][197][198], CCD camerabased systems [199][200][201][202][203][204][205][206] and tomographic system for sub-surface PS distribution measurements only [START_REF] Johansson | Realtime light dosimetry software tools for interstitial photodynamic therapy of the human prostate[END_REF]105,[207][208][209] One of the most efficient approaches to predict PDT damages remains the measurement of singlet oxygen concentration's level [182,193,210]. During a PDT treatment, a luminescence of the singlet oxygen in the near infrared (approximately 1270 nm) appears. Nonetheless, a direct in vivo measurement of this signal remains challenging due to its weak amplitude and the very short lifetime of the molecule (approximately 50 ns). Using photomultiplier tubes [211][212][213] or other technologies [214][215][216], the signal is integrated during a part or the entire treatment time in order to collect enough photons to obtain a relevant metric related to PDT dose. Although this metric correlated closely damages due to PDT treatment, a practical system providing robust results for clinical routine seems no to be available yet.

Additionally, the oxygen partial pressure can also be used as a surrogate metric to correlate the PDT damages. Variations of the partial pressure of oxygen in the blood have an impact on the oxygenation of tissues where oxygen is required for PDT treatment [217][218][219][220].

The study of Li et al. demonstrates a correlation between singlet oxygen production and oxygen partial pressure during irradiation. A real-time multi-parameter detection system integrating PS florescence, luminescence of singlet oxygen, oxygen partial pressure should provide robust data to implement a complete dosimetric assessment and to monitor a PDT treatment efficiently.

Chapter conclusion

PDT dosimetry is a multifactor estimation to compute impact on tissues. Predictive dosimetry includes computing of light, of PS and oxygen distribution in biological tissues.

Although it is highly time-consuming, Monte-Carlo remains the reference method to compute the propagation of light in heterogeneous media. Monitoring of the PDT treatment gathers several features (PS fluorescence, detection of singlet oxygen) to improve the "PDT dose" metric. Finally, a collaborative effort between contributions of biologists, chemists, physicists, and engineers would lead to define a comprehensive PDT dose [221]. In this thesis, a focus was achieved on the light distribution modeling using the Monte-Carlo method. As detailed in chapter II, dosimetry requires an accurate estimation of light dose deposited in the matter. Accordingly, Monte-Carlo model that simulates the photons' interactions in the matter during their entire "life" is the prime candidate to be used in PDT dosimetry. After introducing the theory of the Monte-Carlo model, its implementation in the context of PDT is presented. Additionally, since the Monte-Carlo modelling is highly timeconsuming, acceleration technics are described to make feasible its integration in a treatment planning system dedicated to PDT in clinical routine. Finally, a comprehensive experimental design is proposed to validate the implementation of the Monte-Carlo model achieved during this thesis.

I. Monte-Carlo model

The Monte-Carlo model simulates the propagation of photons through matter according to physics laws. Photons are generated at a random position according to the light source location. Space is discretized according to matrix or grid of voxels and each photon pathway through different media and its progressive absorption are simulated at each voxel location. Thus, by generating a very large number of photons, a light dose (i.e., fluence [J/cm²]) is obtained using the recording of all photons' absorption inside each voxel of the grid.

Theory of the Monte-Carlo model

The Monte-Carlo algorithm describes the history of a photon from its "birth" or launch until its death. As an illustration, in figure 3.1, a photon is launched with an initial state and is propagated through different media. Several scattering and absorption events occur until the photon death. Cartesian coordinates), a direction (ux, uy, uz in Cartesian coordinates) and a weight w (dimensionless). This weight, initially set to 1.0, could be seen as a "health bar" where each collision (i.e., absorption event) decreases this gauge until death (i.e., w = 0.0). To simplify the comprehension of the model, the progression of only one photon will be described in this paragraph. However, each step of this photon's life is common to the others. Only the random characteristic leads to different pathways for each photon generated during the simulation.

𝑝ℎ𝑜𝑡𝑜𝑛 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ ≝ ( ( 𝑥 𝑦 𝑧 ) ( 𝑢 𝑥 𝑢 𝑦 𝑢 𝑧 ) 𝑤 ) (3.1)
All different types of light sources can be modeled using these features (collimated uniform or Gaussian beam, focused or not, isotropic sources). Spherical or cylindrical diffusers are mainly used in interstitial or intraoperative PDT. These light sources are isotropic, i.e., photons are distributed uniformly around the point source. This non-preferential direction of photon direction is initially calculated using eq. 3.2:

𝑢 𝑥 = 𝑠𝑖𝑛 𝜃. 𝑐𝑜𝑠 𝜑 𝑢 𝑦 = 𝑐𝑜𝑠 𝜃. 𝑠𝑖𝑛 𝜑 𝑢 𝑧 = 𝑐𝑜𝑠 𝜃 (3.2)
where 𝜃 is the angle between trajectory vector and the z axis and 𝜑 is the angle between trajectory vector and the x axis.

To obtain an isotropic uniform distribution of photons, 𝜃 and 𝜑 are computed using a random number (RND):

𝑐𝑜𝑠 𝜃 = 2𝑅𝑁𝐷 -1 𝑠𝑖𝑛 𝜃 = √(1 -𝑐𝑜𝑠 2 𝜃) 𝜑 = 2𝜋. 𝑅𝑁𝐷 if 𝜑 < π 𝑠𝑖𝑛 𝜑 = √(1 -𝑐𝑜𝑠 2 𝜑) else 𝑠𝑖𝑛 𝜑 = -√(1 -𝑐𝑜𝑠 2 𝜑) (3.3)
Once the initial position and direction are computed for each photon, their propagation begins. According to the optical properties of the medium, photons cross a given distance called step size sphoton [m]. This step size denotes the displacement of photons through matter at each iteration of the algorithm (i.e., each step of the "photon's life") and is deduced from absorption and scattering coefficients through the probability p(s) (equation 3.4):

𝑝(𝑠) = 𝑒 -𝜇 𝑡 .𝑠 𝜇 𝑡
such as the probability density function is

𝐹(𝑠) = ∫ 𝑒 -𝜇 𝑡 .𝑠 𝜇 𝑡 ∞ 0 𝑑𝑠 = 1 (3.4)
where the total attenuation coefficient µt is the sum of the absorption and scattering coefficients.

Thus, for a given photon and its step size sphoton, the probability distribution function F becomes:

𝐹(𝑠 𝑝ℎ𝑜𝑡𝑜𝑛 ) = ∫ 𝑒 -𝜇 𝑡 .𝑠 𝜇 𝑡 𝑠 𝑝ℎ𝑜𝑡𝑜𝑛 0 𝑑𝑠 = 1 -𝑒 -𝜇 𝑡 .𝑠 𝑝ℎ𝑜𝑡𝑜𝑛 (3.5)
sphoton, describing the iterative photon pathway, is finally obtained by a random sampling of the distribution function F (equation 3.6):

𝐹(𝑠 𝑝ℎ𝑜𝑡𝑜𝑛 ) = 1 -𝑒 -𝜇 𝑡 .𝑠 𝑝ℎ𝑜𝑡𝑜𝑛 = 𝑅𝑁𝐷 𝑠 𝑝ℎ𝑜𝑡𝑜𝑛 = -ln (1 -𝑅𝑁𝐷) 𝜇 𝑡 = -ln (𝑅𝑁𝐷) 𝜇 𝑡 (3.6) Once this first step length is computed, the photon begins its propagation through surrounding matter. Future interactions with media, defined by their optical properties, is iteratively computed from the Hop-Drop-Spin process.

b. Hop-Drop-Spin model description

The Hop-Drop-Spin sequence is the approach used to iteratively compute state and position of the photon along its trajectory. "Hop" designates the displacement of the photon, "Drop" describes to the absorption event and "Spin" defines the change of direction.

i. Hop

During a movement at a given step of its pathway, the photon may change of medium even if its step size is low (only several micrometers). Once the step size is computed, the algorithm evaluates if the photon is changing of medium or not.

If not, a standard photon movement called "Hop" is achieved. The photon moves according to the step size 𝑠 𝑖 , the direction cosines 𝑢𝑥 𝑖 , 𝑢𝑦 𝑖 , 𝑢𝑧 𝑖 and the coordinates xi, yi, zi:

𝑥 𝑖+1 = 𝑥 𝑖 + 𝑠 𝑖 . 𝑢 𝑥,𝑖 𝑦 𝑖+1 = 𝑦 𝑖 + 𝑠 𝑖 . 𝑢 𝑦,𝑖 𝑧 𝑖+1 = 𝑧 𝑖 + 𝑠 𝑖 . 𝑢 𝑧,𝑖 (3.7)
where 𝑖 ∈ [0; 𝑛], n represent the total number of the photon's displacements during its entire history.

Conversely, if the photon moves towards a medium with different optical coefficient, the photon is placed at the boundary of the first medium using the current direction cosines 𝑢 𝑥,𝑖 , 𝑢 𝑦,𝑖 , 𝑢 𝑧,𝑖 and current coordinates xi, yi, zi:

𝑥 𝑖+1 = 𝑥 𝑖 + (𝑠 𝑖 -𝑠 𝑏𝑜𝑢𝑛𝑑𝑎𝑟𝑦 ). 𝑢 𝑥,𝑖 𝑦 𝑖+1 = 𝑦 𝑖 + (𝑠 𝑖 -𝑠 𝑏𝑜𝑢𝑛𝑑𝑎𝑟𝑦 ). 𝑢 𝑦,𝑖 𝑧 𝑖+1 = 𝑧 𝑖 + (𝑠 𝑖 -𝑠 𝑏𝑜𝑢𝑛𝑑𝑎𝑟𝑦 ). 𝑢 𝑧,𝑖 (3.8)
where 𝑠 𝑏𝑜𝑢𝑛𝑑𝑎𝑟𝑦 is the step length between the current position and the boundary with the neighboring medium.

If the neighboring medium owns the same refractive index, no additional interaction can happen. However, when the refractive medium of the current medium 𝑛 1 is different from the surrounding medium 𝑛 2 , the reflectance R must be evaluated to compute probability that the photon could be reflected at the interface or transmitted into the neighboring medium. This reflectance is computed using the Fresnel reflectance equation (equation 3.9):

𝑅 = (sin 𝜃 1 cos 𝜃 2 -cos 𝜃 1 sin 𝜃 2 ) 2 2 × (cos 𝜃 1 cos 𝜃 2 + sin 𝜃 1 sin 𝜃 2 ) 2 + (cos 𝜃 1 cos 𝜃 2 -sin 𝜃 1 sin 𝜃 2 ) 2 (sin 𝜃 1 cos 𝜃 2 -cos 𝜃 1 sin 𝜃 2 ) 2 . (sin 𝜃 1 cos 𝜃 2 + cos 𝜃 1 sin 𝜃 2 ) 2 (3.9)
where 𝜃 1 is the photon incident angle, 𝜃 2 is the photon transmission angle through the second medium computed using the Snell-Descartes law and according to the refractive indices of the first and second media (equation 2.14).

Once the reflectance is computed, a random number between [0,1] is generated and compared to the reflectance value. If the random value is strictly superior to the reflectance value, the photon is transmitted to the second medium with the following direction:

𝑢 𝑥,𝑖+1 = 𝑢 𝑥,𝑖 . 𝑛 1 𝑛 2 (3.10) 𝑢 𝑦,𝑖+1 = 𝑢 𝑦,𝑖 . 𝑛 1 𝑛 2 if cos 𝜃 1 ≈ 1 (normal incidence) 𝑢 𝑧,𝑖+1 = cos 𝜃 1 else if cos 𝜃 1 ≈ 0 (raking incidence) 𝑢 𝑧,𝑖+1 = 0 else 𝑢 𝑧,𝑖+1 = cos 𝜃 2
If the random value is inferior to the reflectance value, the photon is reflected to the interface between the two media with the following direction cosines:

𝑢 𝑥,𝑖+1 = 𝑢 𝑥,𝑖 𝑢 𝑦,𝑖+1 = 𝑢 𝑦,𝑖
𝑢 𝑧,𝑖+1 = -𝑢 𝑧,𝑖 (3.11) ii. Drop

Once the "Hop" movement completed, the "Drop" (i.e., absorption) event occurs.

After moving to its new position, the photon interacts with the tissue. A small proportion of its weight w [dimensionless] is decreased according to the size step 𝑠 𝑖 :

𝑤 𝑖+1 = 𝑤 𝑖 . 𝑒 -𝑠 𝑖 .𝜇 𝑎 (3.12)
The absorbed fraction of the photon by the medium has to be recorded in order to obtain the total fluence deposit by all photons generated. To accumulate these absorbed fractions, the entire space is discretized into voxels which generates a matrix called Absorption.

The absorption is then stored in the Absorption matrix at position l, m, n, which corresponds to the location where the absorption occurs:

𝐴𝑏𝑠𝑜𝑟𝑝𝑡𝑖𝑜𝑛[𝑙, 𝑚, 𝑛] 𝑖 = 𝐴𝑏𝑠𝑜𝑟𝑝𝑡𝑖𝑜𝑛[𝑙, 𝑚, 𝑛] 𝑖-1 + 𝑤 𝑖 . (1 -𝑒 -𝑠 𝑖 .𝜇 𝑎 ) (3.13)
iii. Spin

Once the "Drop" event completed, the photon is scattered into a new direction according to the optical properties of the medium, which corresponds to the "Spin" process.

These new direction cosines are computed using the Henyey-Greenstein equation seen previously (equation 2.7). The function can be normalized over 4π sr (equation 3.14):

∫ {∫ 𝑝(θ). 𝑠𝑖𝑛(θ). 𝑑θ π 0 } 2π 0 𝑑𝜙 = 1 (3.14)
Thus, the Henyey-Greenstein function becomes:

𝑝(θ) = 1 2 . 1 -𝑔² (1 + 𝑔 2 -2𝑔. cos (θ)) 3 2 ⁄ (3.15)
The new direction is obtained by inserting a random number RND in the previous equation (equation 3.15). New angles θ and 𝜑 provide from the sampling of the Henyey-Greenstein function:

cos θ = 1 + 𝑔 2 -( 1 -𝑔 2 1 -𝑔 + 2𝑔. 𝑅𝑁𝐷 ) 2 2𝑔 𝜑 = 2𝜋. 𝑅𝑁𝐷 (

c. Termination of the history

The photon keeps propagating using the hop-drop-spin sequence with a decrease of its weight w at each interaction. The roulette method is one of the most common used to finish the photon propagation. An arbitrary threshold value is defined (10 -4 typically). When the photon's weight reaches a value lower than this threshold, the roulette algorithm is applied. A random number RND is generated and compared to a chance value (typically 0.1). If the RND value is strictly superior to the chance value, the photon terminates its propagation.

Conversely, if the RND value is inferior to the chance value, the photon carries on with a weight value increased by the inverse of the chance value (typically 10 times the last weight value).

𝑖𝑓(𝑤 𝑖 < 𝑡ℎ𝑟𝑒𝑠ℎ𝑜𝑙𝑑) 𝑖𝑓(𝑅𝑁𝐷 ≤ 𝑐ℎ𝑎𝑛𝑐𝑒) 𝑤 𝑖+1 = 𝑤 𝑖 𝑐ℎ𝑎𝑛𝑐𝑒 𝑒𝑙𝑠𝑒 𝑤 𝑖+1 = 0 (3.19)
Once the photon history is terminated, another photon is launched. The Monte Carlo simulation is repeated until the expected total number of photon Nphotons is reached. For example, to compute 10 6 photons in a homogeneous medium composed with the absorption µa = 0.2 cm -1 and scattering µs = 160 cm -1 coefficients, which correspond to mean brain optical properties, (a mean tstep with such coefficients is approximately 700 ns2 ) and with a threshold fixed to 10 -4 , the required mean time is about 65 minutes.

Finally, all these previously described steps are summarized in figure 3.2. A photon is launched and its step size is computed according to optical coefficients of the medium. Then, the hop-drop-spin sequence is computed to move through matter until the photon's death.

This process is then repeated to reach a large number of generated photons. Graaff et al. in 1993 [226].

ii. Perturbation methods

Based on a similar idea, the perturbation method requests previous Monte-Carlo simulations to be achieved [227][228][229][230][231][232][233][234][235][236][237][238][239] to consider resolution in a heterogeneous medium.

Several features deduced from previous Monte-Carlo simulation in homogeneous media are stored all along the path length, the exit weight and the number of interactions encountered during all photons histories. After collecting these data, inhomogeneities are introduced into the structure of the second Monte-Carlo simulation and impact the photon weight [227,230,233,240].

iii. Hybrid methods

Another approach to accelerate Monte-Carlo simulation is to introduce analytical solutions in determined areas in order to avoid a full Monte-Carlo simulation on the whole space.

i. CPU parallelism

The first published attempt to parallelize algorithm was Kirkby et al. in 1997 [255].

The authors proposed to build a network between several computers and parallelize Monte-Carlo simulation on these different CPU. Once 1000 photons were generated on all the computers, results were stored in a file on a disk that was accessible to all connected computers in the network. A maximum of 24 computers were declared because each of them

II. GPU Monte-Carlo implementation

Monte-Carlo method is considered as the reference in the dosimetry computation for both ionizing and non-ionizing radiation. Because this technic is highly time-consuming in conventional central processing unit (CPU) implementation, several improvements have been developed to speed-up simulations. Among them, the use of GPU in order to parallelize algorithm provides a cost-efficient solution to accelerate Monte-Carlo method. In this paragraph, a parallel implementation of Monte-Carlo using GPU technology is described.

Monte-Carlo implementation

a. History

One of the most important methods for the computer performances improvement has been, for four decades, to increase the speed at which the processor's clock operated.

Increasing this clock speed resulted in also increasing the number of operations performed per second by the CPU. The multiplication of CPU in supercomputer has been the starting point of the new point of view in the escalation of computer power. Since 2005, instead of stretching the limits of processor's clock speed, the use of multi-core CPU (2 to 16 cores) has been allowing to strongly increase performances by parallelizing tasks. In 2006, the manufacturer NVIDIA unveiled the first GPU which was not dedicated to produce a feed of output images to a display but to enable general-purpose computation. Thus, this new electronic component could be seen as a multi-core CPU owning thousands of tiny processing units. Currently, general-purpose computing on graphics processing units (GPGPU) is on the rise in image processing, fluid dynamics and bioinformatics.

b. Introduction to parallel world

The main goal of parallelized implementation is to create an algorithm capable to be executed simultaneously in different processors in order to obtain the output in less time. To exemplify this concept, matrix calculation can be easily integrated in a parallel computing.

Let us consider the sum of the two 𝑚 × 𝑛 matrices A and B: To be executed in parallel, a computer code must be cut into different independent processes.

[ 𝑎 1,1 ⋯ 𝑎 1,𝑛 ⋮ ⋱ ⋮ 𝑎 𝑚,1 ⋯ 𝑎 𝑚,𝑛 ] + [ 𝑏 1,1 ⋯ 𝑏 1,𝑛 ⋮ ⋱ ⋮ 𝑏 𝑚,1 ⋯ 𝑏 𝑚,𝑛 ] = [ 𝑎 1,1 + 𝑏 1,1 ⋯ 𝑎 1,𝑛 + 𝑏 1,𝑛 ⋮ ⋱ ⋮ 𝑎 𝑚,1 + 𝑏 𝑚,1 ⋯ 𝑎 𝑚,𝑛 + 𝑏 𝑚,𝑛 ] (3.22
Several criteria have been defined to evaluate performances. Total computing time 𝑇 is the sum of the sequential computing time 𝑇 𝑆 and parallel computing time 𝑇 𝑃 : 𝑇 = 𝑇 𝑆 + 𝑇 𝑃 (3.23) According to the Amdhal's law [276,277], the minimum computing time 𝑇 𝑛𝑝 on a computer owning np processing units is equal to:

𝑇 𝑛𝑝 = 𝑇 𝑆 + 𝑇 𝑃 𝑛𝑝 (3.24)
The acceleration 𝐴 𝑛𝑝 due to the parallelism is then defined as:

𝐴 𝑛𝑝 = 𝑇 𝑇 𝑛𝑝 = 1 (1 -𝛼) + 𝛼 𝑛𝑝 (3.25)
Where 𝛼 is the proportion of execution time that benefits from parallelized implementation (e.g., according to a given code implemented, if 60% of its execution time could be subject of a speedup, 𝛼 is equal to 0.6).

The parallelization efficiency 𝐸 𝑛𝑝 is the ratio of the acceleration 𝐴 𝑛𝑝 and the number of processing units np:

𝐸 𝑛𝑝 = 𝐴 𝑛𝑝 𝑛𝑝 = 1 (1 -𝛼). 𝑛𝑝 + 𝛼 (3.26)
As an example, if 𝛼= 0.99, which means that 99% of the entire programming code could be parallelized on a computer owning np = 99 processing units, the acceleration is equal to 50 and parallelization efficiency is equal to 0.51. Thus, the one percent of the programming code that cannot be parallelized takes the half time of the entire computing time.

Two rules must be validated in order to know whether a programming code can be parallelized or not. Firstly, the result must not be altered if a random permutation occurs between processing unit execution. If an inversion is made provoking the execution of the code in the unit execution np instead of the unit execution np-1, the output must not be affected. Secondly, two codes in execution must be totally independent. The output of the unit execution np must be not the entry of the unit execution np-1 (data dependence) and the output of the two units execution np-1 and np must not change the same variable (output dependence).

c. Architecture

The particularity of the GPGPU programming lies in the architecture of the device used. The term "Host" refers to the CPU and its memory; the "Device" to the GPU. Only the parallelized code is equally distributed into several unities called "kernel" as referred previously as the "tiny processing units" and located in the device (see figure 3.4). Each kernel instance is called "thread". Several kernels constitute a "block" and several blocks constitute a "grid". In some cases, it can be useful to define several grids. It is commonly used to illustrate this interlocking in 2D or 3D. However, no determined organization of kernels is defined in the device. This compartmentalization of the GPU memory is defined according to hierarchy of electronic components embedded in the printed circuit board of the graphic processor. While the CPU memory can be considered as unique, the GPU memory is divided into six different memories as follows (see figure 3.5):

-Global memory: this is the largest memory available in the device. The main asset is that stored data are visible from all kernels in the device and can store large amount of data. Each thread can read and write in this memory but the access is very slow.

-Constant memory: this one holds only data that will remain constant during the execution of the code. Thus, its authorization is in read-only, which allows a more rapid access to constant data.

-Texture memory: similar to the constant one, it can be used when memory access exhibits a strong dependence with local neighbor threads. When all read accesses in a warp (merging of 32 threads) are physically adjacent, texture memory can decrease memory traffic and increase performances in comparison to the global memory.

These next three memory types are visible to the host which can read and write.

-Local memory: it is contained into each block and is visible only to one thread.

Threads have read and write access to this memory.

-Register: this one owns the same specifications as the local memory. The space available is smaller than the local memory but, because registers reside on the GPU chip, the access to data is very fast.

-Shared memory: this memory resides also on the GPU chip but is visible to all threads contained in a block in a fast read and write access.

Using all these six different memories efficiently, program performances can be optimized and the parallelization efficiency improved. The matrix corresponding to the segmentation output (i.e., biological medium of each voxel) is also stored in constant memory. This matrix is composed of Nx.Ny.Nz indices. Each index refers to a lookup table in which optical properties of each medium type are stored.

Each thread generated on the GPU runs an entire photon history (from its birth with initial properties to its death by the roulette process) (see figure 3.6). Thus, depending on hardware capacities, several photon histories are simultaneously launched. Propagation process is the same as presented in chapter II, namely:

-A photon is first created at a random position between the two extremities (xin, yin, zin), This PRNG provides a lighter implementation than the classical Mersenne Twister algorithm and is optimized to GPGPU requiring large state generators. Moreover, it has been evaluated using the strict statistical test "TestU01" provided by L'ecuyer et al. [283] and results reported pretty good quality of the outputs. Although, the initial Tiny Mersenne Twister own a period of generated sequences of 2 19937 -1 iterations, the tiny version owns one of 2 127 -1. Besides its correct performances, this particular PRNG has been selected and implemented in the Monte-Carlo mostly because it is available on OpenCL and CUDA platforms.

The PRNG is seeded in each thread using the time corresponding to the program's launch and the index of the thread in order to obtain a different seed in every thread.

b. Atomic operation

The atomicity property is the second critical point in GPU Monte-Carlo. This term designates an operation (in most cases addition, subtraction, multiplication or division)

performed by two or more threads on the same allocated memory location and at the same time. The probability of such an addressing conflict is important due to the numerous threads created in a GPU program. In the global memory, which is visible for every thread, several threads can expect to write in the same allocated memory location. When a large number of threads read and write on the same global memory location concurrently, which is the case in the "drop" event typically, part of concurrent threads will not succeed to write in the global memory. If this property is not taken into account, absorption can be underestimated.

To circumvent this issue, atomic operations have been created to allow this accumulation on the same memory location without any loss and to ensure the integrity of results. For example, if a thread reads a data at a given address in global or shared memory, modifies the contained value and writes back the result to the same address, no other thread can access this address during this operation. Thus, this modification is guaranteed to be performed without competitive issue from other threads. However, these operations are computationally expensive and slow down the acceleration due to waiting threads during an atomic process.

c. Synchronization

Because photons propagation in Monte-Carlo method is governed by random properties, a photon propagated in a thread i and another photon propagated in the thread i+1

will not terminate their history at the same time. To counterbalance this probable time interval, the device can synchronize all the threads that are contained in the same block. When the thread i reaches the synchronize command, it stops its computing until all threads contained in its block also reaches this command. Thus, this barrier allows to keep a homogeneous progression of the threads during the simulation.

d. Memories optimization

The last step to improve performances and rapidity of the execution is the optimization of the memory's allocation. Although the constant memory is relatively small (i.e., 65536

bytes, approximately 8 Mo), it remains sufficiently large to store the 12 input variables ( 48octets). Conversely, matrix containing indices of tissues owns the same dimensions as MRI standard volume, i.e., 256x256x150 voxels. This matrix is approximately 40Mo and cannot be sent into the constant memory. Also, because absorption matrix will be constantly modified during the execution and owns the same dimensions as the tissue matrix, these two matrixes must be stored into the global memory. Texture memory cannot be used because there is no spatial dependence between threads.

Regarding the number of grids, blocks and threads per block, no law is currently available to compute the optimal configuration. Several tests must be performed to evaluate the best configuration. Moreover, this arrangement is only optimal for a given GPU. That is why a list of recommended GPU is often provided with particular GPU program.

III. Model validation

Once the GPU Monte-Carlo implemented, a comparison with experimental measurements is expected to validate results. This paragraph details the methodology used to assess the results of Monte Carlo algorithm from experimental measurements. Additionally, optimization of the performances to produce fast Monte-Carlo simulations (such as memories optimization and programming language) is highlighted.

Phantom creation

A phantom of known optical properties is firstly created. Several types of phantom are used to mimic precise optical properties, such as optical polymer gels [284][285][286], silicone [287,288],

Al2O3 powders [289,290] or 3D printing [291,292] However, the most popular and easier way to build such a phantom is to mix Intralipid liquid, India ink and water [194,[293][294][295][296]. This type of phantom is usually applied as a biological Values of absorption µa = 0.02 mm -1 and reduced scattering µ's = 2 mm -1 correspond to mean optical properties of normal human brain tissues with infiltrated glioma cells [START_REF] Stepp | ALA and Malignant Glioma: Fluorescence-Guided Resection and Photodynamic Treatment[END_REF]117,297,298]. These values will be considered as references in this study.

Thus, for Ԑa,ink = 384 ± 8 mm -1 , Ԑ's,il = 26.3 ± 0.4 mm -1 , µa = 0.02 mm -1 , µ's = 2 mm -1 and mtot = 1000 mg, masses of India ink, Intralipid 20% and water are respectively mink = 0.052g, mil = 76.046g and mwater = 923.902g. a calibration factor CFmedium dependent of the surrounding medium must be calculated to enable this conversion:

Monte-Carlo validation by comparison with experimental measurements: materiel and methods

A

ϕ medium = CF 𝑚𝑒𝑑𝑖𝑢𝑚 . P measured (3.28) In their two papers, Marijnissen and Star provide a complete methodology to obtain this calibration factor CFmedium [299, 300] which is recommended by the American Association of Physicists in Medicine (AAPM) [301,302]. This factor is the combination of a calibration factor determined firstly in the air CFair and four other corrective coefficients.

CF 𝑚𝑒𝑑𝑖𝑢𝑚 = CF 𝑎𝑖𝑟 . F 𝑛 . F 𝑝 . F 𝑏 . F 𝑖 (3.29) where F 𝑛 is a correction factor enabling to compensate the loss of light due to interface between the sensor and the surrounding medium using their respective refractive indices

[dimensionless]. F 𝑝 is a perturbation correction [dimensionless].
F 𝑏 is a correction of the collection photons loss due to the surface occupied by the fiber on the detection sphere [dimensionless].

F 𝑖 is a correction of the sensor non-homogeneity response [dimensionless].

Combining these empirical and theoretical factors, a calibration factor CF 𝑚𝑒𝑑𝑖𝑎 can be estimated for a given biological medium and for a given isotropic probe in order to obtain fluence rate value ϕmedium [W.cm -2 ] in this medium:

ϕ medium = P measured . CF 𝑎𝑖𝑟 . F 𝑛 . F 𝑝 . F 𝑏 . F 𝑖 (3.30) First, a calibration factor in air CF 𝑎𝑖𝑟 must be determined (see figure 3.7). A bare fiber was fixed on an experimental optical bench and connected to a 630nm laser source. In the other side of the optical bench, two detectors were placed in front to the output of the bare fiber. The first sensor is an isotropic probe (IP85, Medlight, Ecublens, Switzerland) which will be used in future power measurements. The sensor was fixed on a rotating support enabling to vary the angle θ between incident beam and tip axe. The second is a flat sensor (818-SL, Newport, Irvine, CA, USA) whose size of the detection cell is 1cm². Direct measurement with this last sensor enables to obtain an estimation of the irradiance [W.cm -2 ].

Thus, measuring the power collected at the same location with both the isotropic probe (a mean is performed according to the irradiation angle θ) P mean_isotropic and the flat sensor, a calibration factor in air CF 𝑎𝑖𝑟 can be estimated: , f = 0.3516 and h = 0.1933 Using these coefficients and for a wavelength of 635 nm, the corresponding anisotropy used was approximately 0.733.

In order to obtain a reduced scattering coefficient µ 𝑠 ′ of 20 cm -1 , a scattering coefficient was deduced using the anisotropy value previously computed:

µ 𝑠 = µ 𝑠 ′ (1 -𝑔) (3.33)
Thus, a homogeneous medium has been implemented with the following optical properties:

-Absorption coefficient: µa = 0. Ecublens, Switzerland and Biolitec GmbH, Jena, Germany) represents the standard optical fiber dimensions used in interstitial PDT.

Monte-Carlo validation by comparison with experimental measurements: results

a. Calibration factor

The calibration factor CF 𝑚𝑒𝑑𝑖𝑢𝑚 was estimated using the following values of correction factors F 𝑛 , F 𝑏 , F 𝑝 and F 𝑖 (see table 3.2 and figure 3.11). This value is only valid for the isotropic probe and the phantom that will be used in the future experimentations.

Correction factors (standard deviation)

Fb (dimensionless) Power measurements were normalized according to each maximum power.

Using the calibration factor previously described, each power measurements have been converted into fluence rate values. Mean relative deviations were obtained according to distances to the fiber tip (see figure 3.13). By computing the mean of all means relative deviations for each cylindrical diffusers and independently from the distance to the light source, a total mean value of all experimentations of approximately 14.97 % ± 6.16 was computed (see table 3.3).

Manufacturer and model of the fiber Mean relative deviations (standard deviation)

Biolitec, CD 603-20

Performances

All the following tests were realized with a NVIDIA GPU Quadro K620. This GPU has the following features (see table 3 -Absorption coefficient: µa = 0.1 cm -1 -Scattering coefficient: µs = 100 cm -1 -Anisotropy coefficient: g = 0.8 102400 photons were generated to simulate a 40 mm cylindrical diffuser placed at the center of a 40x40x80 mm homogeneous slab constituted of the medium previously described.

First, the number of blocks has been varied from 1 to 1024. The number of threads per block was fixed to 512 during all simulations. Mean computing time of the GPU part only was recorded. Results are presented in the following table 3 Strong dependence of computing time is observed according to the number of blocks employed. In the case of the GPU previously described and used here, the number of 32 blocks appears as the best solution to improve GPU performances.

Once the number of 32 blocks was fixed, the number of threads per block has been varied from 8 to 512. The same simulation conditions were implemented.

b. Fluctuations

The minimal photons' number that must be generated to obtain an accurate result (less than ten percent of error with a very large photon number) was evaluated. The number of photons was varied from 10 5 to 10 7 using the same properties as defined previously. A homogeneous medium was implemented with the following optical properties:

-Absorption coefficient: µa = 0.1 cm -1 -Scattering coefficient: µs = 100 cm -1 -Anisotropy coefficient: g = 0.8

A 30 mm cylindrical diffuser was implemented at the center of this slab. Fluence rate values ϕ were extracted along the cross section of the cylindrical diffuser's center. Relative error was computed at several distances d from the center of the source. The reference simulation used 10 8 photons. Five different total numbers of photons (10 4 , 10 5 , 10 6 , 10 OpenCL platform using the same computing parameters.

Although the OpenCL platform enables all GPU as device, development and optimization is much more challenging than CUDA. Moreover, using the same algorithm with the same computing parameters, OpenCL reveals to be more than seven times slower than the CUDA simulation (see table 3.8).

Discussion and chapter conclusion

Several different approaches enable to compute the propagation of light. Among them, Monte-Carlo remains the reference method. Despite its high accuracy, this method remains so far little used mainly because of the long computing time required to perform a whole simulation. With recent high computing capacities and technologies enabling the parallel computing, Monte-Carlo method can now reach a computing time acceptable for clinical routine.

This chapter introduces a programming code fully dedicated to the GPU modeling of cylindrical diffusers used in PDT is described. To evaluate the output of the algorithm, i.e., the fluence rate emitted from light sources, a phantom has been created to mimic particular optical properties. Experimental results were compared to Monte-Carlo simulations and a mean relative deviation of approximately 14.97% enabled to validate the accuracy of the Monte-Carlo model.

Although the value of the total mean could appear as high, it fits with the recommendations of the AAPM [301,302] and enables to validate the accuracy of the Monte-Carlo model proposed:

"the potential variations in the many factors that affect PDT can easily exceed a minor 10% to 20% error in light dosimetry." [302] Several biases could be noticed to explain these mean relative deviations:

-The precision of the weighting of every compound to create the optical phantom is critical. Particularly, the mass of India ink injected into the phantom is low but has a strong impact on absorption coefficient. A slight measurement error during the weighting of the India ink can provoke a shift in the phantom absorption.

-The precision of the isotropic probe location in the phantom according to the cylindrical diffusers position was evaluated manually. Once the liquid optical phantom was poured inside the black box, the isotropic probe could not be seen. It was then impossible to check the position of the bulb during the experiment.

-The efficiency of the cylindrical diffusers used has been evaluated using three different laser systems: Ceralas (Biolitec, Jena, Germany), DIOMED 630 PDT (DIOMED, Cambridge, United Kingdom) and ML7710-630-6K (Modulight, Tampere, Finland). The mean of the three efficiencies measured was used to weight the Monte-Carlo simulations.

-The inhomogeneity of illumination profiles of the cylindrical diffusers was measured on each diffuser tip and shows a strong diversity of illumination profiles. Ideally, this heterogeneity should be measured and considered during planning therapy in order to reach a more relevant dosimetry.

-The optical coefficients used in the Monte-Carlo simulations to simulate the phantom are based on forecasting that the phantom is composed of these same coefficients.

-The non-realistic modeling of a cylindrical diffuser in the Monte-Carlo simulations impacts the results. The non-diffusing part and the numerical aperture of each cylindrical diffuser should be considered to obtain a more realistic modeling. These characteristics influence the light distribution.

In clinical context, the strong optical heterogeneity of brain tissues, and particularly in the GBM, induces multiple perturbations for light propagation. Also, optical coefficients provided by different studies [305,306] For non-resectable GBM, newly diagnosed or relapse, interstitial photodynamic therapy (iPDT) may be preferred as alternative treatment modality. This surgical procedure aims for inserting light sources directly inside the tumor and requires the use of a treatment planning system (TPS) before entering the operating theater.

A proof of concept of a TPS dedicated to the stereotactic iPDT treatment of GBM is introduced in this chapter. The two key points in the procedure planning are the stereotactic localization and the dosimetry. These two sections are evaluated and the respective methodologies are described.

Since this first version is only a proof of concept, forthcoming developments are then proposed, particularly the integration of the GPU Monte-Carlo method presented in chapter III.

I. A treatment planning system

For unresectable GBM cases, iPDT may be preferred [19,117,118,297,308]. Since this procedure requires to plan the positioning of light sources to compute the light exposure duration, a TPS dedicated to this modality is highly expected.

Interstitial PDT procedure

According to the literature, a minority [START_REF] Paquette | Prognostic value of health-related quality of life for death risk stratification in patients with unresectable glioblastoma[END_REF] to 23%-38% [START_REF] Fazeny-Dorner | Survival and prognostic factors of patients with unresectable glioblastoma multiforme[END_REF] of GBM are nonresectable (newly diagnosed or relapsing). In these cases, iPDT [117,[308][309][310][311]

Proof of concept

A TPS is highly recommended to plan an iPDT procedure and for optimization of the light delivery and maximization of the treated volume. Once this modeling is completed, the illumination duration is computed in order to deposit a given therapeutic fluence value [J.cm After a target is manually delineated, the surgeon or medical physicist is able to simulate different treatment configurations through the placement of cylindrical diffusers inside the target volume. To simulate the fibers insertion, the diffuser length [mm], the power of the source [mW], and the input and output points (i.e., entry point into the skull and ending point into the brain) must be defined. Vascular structures and optic fibers collision is a major concern in this stereotactic procedure. Specific visualization tools are thus required to enable visual inspection of the fiber trajectory and thus limit damages to the vascular tree. The fiber trajectory within the volume may be represented in two ways (see figure 4.3). When the whole trajectory is contained in only one axial image, it is represented by two distinct parts: the diffusing part, symbolized here by a red continuous line, and the non-diffusing part, symbolized here by a cyan continuous segment. Additionally, when the trajectory crosses several planes, the projection of the trajectory on the displayed plane is symbolized here by a discontinuous segment navy blue (see figure 4.3). Thus, it is possible to visualize the progression of the fiber inside the patient's brain. Entry and ending points are directly Leksell frame were included inside at least two sides of the box. Thus, these six rods appear as points in the axial plane of the patient imagery (see figure 4.4). To perform a Leksell localization of the DICOM set, these rods have to be located on all images. Two collections of points in 𝑅 3 were created: the first 𝑃 𝐼𝑀𝐴 represents four points on n slices in the imaging modality (i.e., CT or MRI) coordinates, and the second 𝑃 𝐿𝑒𝑘𝑠𝑒𝑙𝑙 represents the theoretical positions of the same four points in the Leksell frame coordinates. A matrix A allows the conversion of 𝑃 𝐼𝑀𝐴 in 𝑃 𝐿𝑒𝑘𝑠𝑒𝑙𝑙 :

𝑃 𝐿𝑒𝑘𝑠𝑒𝑙𝑙 = 𝐴 . 𝑃 𝐼𝑀𝐴 (4.1)
For a given slice, the four points enclosed in red circles (see figure 4.4) form a rectangle 190 x120 mm in length and the intersection of diagonals is situated on x = y = 100 mm Leksell coordinates. They also own the same x and y Leksell frame coordinates. The two points enclosed in a blue circle (see figure 4.4) indicate the location of the diagonal rod of the Leksell frame. These two points move, respectively, along the lines [1,2] and [3,4] according to the axial imaging modality slice location. Thus, the z Leksell coordinate can be determined PIMA and PLeksell are thus composed of 4n 3D (x,y,z) elements: (4.4) With these two vectors, a 3x3 covariance matrix S is built: 𝑆 = ( 𝐶𝑜𝑣(𝑃 𝐼𝑀𝐴 (: ,1), 𝑃 𝐿𝑒𝑘𝑠𝑒𝑙𝑙 (: ,1)) 𝐶𝑜𝑣(𝑃 𝐼𝑀𝐴 (: ,2), 𝑃 𝐿𝑒𝑘𝑠𝑒𝑙𝑙 (: ,1)) 𝐶𝑜𝑣(𝑃 𝐼𝑀𝐴 (: ,3), 𝑃 𝐿𝑒𝑘𝑠𝑒𝑙𝑙 (: ,1)) 𝐶𝑜𝑣(𝑃 𝐼𝑀𝐴 (: ,1), 𝑃 𝐿𝑒𝑘𝑠𝑒𝑙𝑙 (: ,2)) 𝐶𝑜𝑣(𝑃 𝐼𝑀𝐴 (: ,2), 𝑃 𝐿𝑒𝑘𝑠𝑒𝑙𝑙 (: ,2)) 𝐶𝑜𝑣(𝑃 𝐼𝑀𝐴 (: ,3), 𝑃 𝐿𝑒𝑘𝑠𝑒𝑙𝑙 (: ,2)) 𝐶𝑜𝑣(𝑃 𝐼𝑀𝐴 (: ,1), 𝑃 𝐿𝑒𝑘𝑠𝑒𝑙𝑙 (: ,3)) 𝐶𝑜𝑣(𝑃 𝐼𝑀𝐴 (: ,2), 𝑃 𝐿𝑒𝑘𝑠𝑒𝑙𝑙 (: ,3)) 𝐶𝑜𝑣(𝑃 𝐼𝑀𝐴 (: ,3), 𝑃 𝐿𝑒𝑘𝑠𝑒𝑙𝑙 (: ,3)) ) (4.5) Using the trace TrS of the S matrix, a 4x4 matrix Q is created from these two data sets:

𝑃 𝐼𝑀𝐴 = { 𝑥 𝐼𝑀𝐴,
𝑄 = [ 𝑇𝑟 𝑆 𝑆 1,2 -𝑆 2,1 𝑆 2,0 -𝑆 0,2 𝑆 0,1 -𝑆 1,0 𝑆 1,2 -𝑆 2,1 2 × 𝑆 0,0 -𝑇𝑟 𝑆 𝑆 0,1 + 𝑆 1,0 𝑆 0,2 + 𝑆 2,0 𝑆 2,0 -𝑆 0,2 𝑆 1,0 + 𝑆 0,1 2 × 𝑆 1,1 -𝑇𝑟 𝑆 𝑆 1,2 + 𝑆 2,1 𝑆 0,1 -𝑆 1,0 𝑆 2,0 + 𝑆 0,2 𝑆 2,1 + 𝑆 1,2 2 × 𝑆 2,2 -𝑇𝑟 𝑆 ] (4.6)
Eigen vectors are computed from the Q matrix and enable the building of a quaternion q expressed as q = a + b + c + d where a, b, c, d are normalized real parts of the maximum eigen values from the eigen vector. With this quaternion, a matrix A which enables the conversion of image coordinates into Leksell coordinates is created:

𝐴 = [ 𝑎 2 + 𝑏 2 -𝑐 2 -𝑑² 2𝑏𝑐 -2𝑎𝑑 2𝑎𝑐 + 2𝑏𝑑 𝑡 𝑥 2𝑎𝑑 + 2𝑏𝑐 𝑎 2 -𝑏 2 + 𝑐 2 -𝑑² 2𝑐𝑑 -2𝑎𝑏 𝑡 𝑦 2𝑏𝑑 -2𝑎𝑐 2𝑎𝑏 + 2𝑐𝑑 𝑎 2 -𝑏 2 -𝑐 2 + 𝑑² 𝑡 𝑧 0 0 0 1 ] (4.7)
where: This fluence rate estimation is computed on each voxel of the space. The coefficients used are provided from a previous study of Beck et al. in 2007 [117] and confirmed by a more recent study of Tedford et al. in 2015 [317]. They correspond to a normal brain tissue infiltrated by GBM tumor cells: µa = 0.2 cm -1 and µs' = 20 cm -1 .

𝑡 𝑥 = 𝑃
To evaluate this analytical model, the results were compared to the Monte-Carlo model presented in chapter III.

To compare these models, a common metric (equation 4. Several Monte-Carlo simulations were computed to evaluate the fluence rate distribution from different diffuser lengths (10, 15, 20, 30, 40 and 50 mm). Metric values were extracted each 1 mm from the center of the light source to a distance of 10 mm (figure 4.5). Analytical model was implemented in Matlab software (MathWorks, Natick, USA).

First, the mean relative deviations for each distance d (between 0.75 mm to 10 mm) were computed over six different diffuser lengths. Subsequently, the mean relative deviations for all six different diffusers were computed over each distance d (between 0.75 mm to 10 mm).

d. Sensitivity study

Fluence rates computed from different analytical models were expected to be very sensitive to optical coefficients of the considered media, particularly when the strong heterogeneity of GBM tissues led to a wide range of optical coefficients (see table 4.5). This sensitivity was evaluated with indices called the Sobol indices Si [318,319].

The estimation of the first-order Sobol index of input parameters Xi requires evaluating the ratio between conditional and total variance (equation 4.12).

𝑆 𝑖 = 𝑉 𝑖 V(𝑌) = V(E[Y|𝑋 𝑖 ]) V(𝑌) (4.12)
This index quantifies the sensitivity of the output Y due to the variation of the input Xi.

𝑉 𝑖 = V(E[𝑌|𝑋 𝑖 ]) = E[E[𝑌|𝑋 𝑖 ] 2 ] -E[E[𝑌|𝑋 𝑖 ]] 2 = 𝑈 𝑖 -E[𝑌] 2 (4.14)
A randomization method estimates Ui using typical variance computing using two Ndimension random samples of input parameters 𝑋 ̃(𝑁) (1) and 𝑋 ̃(𝑁) (2) , with Xi constant: 1) , … , 𝑥 𝑘(𝑖-1)

𝑈 ̂𝑖 = 1 𝑁 ∑ 𝑓(𝑥 𝑘1 (
(1)

, 𝑥 𝑘𝑖 (1) , 𝑥 𝑘(𝑖+1) (1) , … , 𝑥 𝑘𝑝 (1) )

𝑓(𝑥 𝑘1 (2) , … , 𝑥 𝑘(𝑖-1)

(2)

, 𝑥 𝑘𝑖 (1) , 𝑥 𝑘(𝑖+1) (2) , … , 𝑥 𝑘𝑝 (2) ) 𝑁 𝑘=1 (4.15) First-order Sobol indices 𝑆 ̂𝑖 can then be estimated as: The same process was applied to estimate the second-order Sobol index; two N-dimension random samples of input parameters 𝑋 ̃(𝑁) (1) and 𝑋 ̃(𝑁) (2) are injected in the f function with Xi and Xj constants: 1) , … , 𝑥 𝑘(𝑖-1)

𝑆 ̂𝑖 = 𝑉 𝑖 V ̂= 𝑈 ̂𝑖 -𝑓 0 ̂2 𝑉 ̂ (
𝑈 ̂𝑖 = 1 𝑁 ∑ 𝑓(𝑥 𝑘1 (
, 𝑥 𝑘𝑖 (1) , 𝑥 𝑘(𝑖+1)

, … , 𝑥 𝑘(𝑗-1)

, 𝑥 𝑘𝑗 (1) , 𝑥 𝑘(𝑗+1)

, … , 𝑥 𝑘𝑝 (1) )

𝑓(𝑥 𝑘1 (2) , … , 𝑥 𝑘(𝑖-1)

, 𝑥 𝑘𝑖 (1) , 𝑥 𝑘(𝑖+1)

, … , 𝑥 𝑘(𝑗-1)

, 𝑥 𝑘𝑗 (1) , 𝑥 𝑘(𝑗+1)

, … , 𝑥 𝑘𝑝 (2) ) 𝑁 𝑘=1 (4.17)

Second-order Sobol indices 𝑆 ̂𝑖𝑗 can then be estimated as:

𝑆 ̂𝑖𝑗 = 𝑉 𝑖𝑗 V ̂= 𝑈 ̂𝑖𝑗 -𝑓 0 ̂2 -𝑉 ̂𝑖 -𝑉 ̂𝑗 𝑉 ̂ (4.18)
As reported by Zhang et al. [319], the parameters and their limits must be carefully defined. In this study, two parameters were included in the sensitivity analysis: absorption and reduced scattering coefficient. Thus, two N-dimension samples (𝑋 ̃(𝑁) (1) and 𝑋 ̃(𝑁) (2) ) of absorption and reduced scattering coefficients were created. Random values of optical coefficients were generated with a range of ± 50% from the mean value presented previous (i.e., µa = 0.02 mm -1 and µs' = 2 mm -1 ): between 0.01 to 0.03 mm -1 for the absorption coefficient and 1 to 3 mm -1

for the reduced scattering coefficient. The Sobol function provided in the "sensitivity" package was applied with a sample size of N = 100,000 and 1000 bootstrap replications. First, second and total-order sensitivity indices were generated and analyzed. The analytical model was implemented with a 40 mm diffuser length. Indeed, the diffuser length can be considered as a constant because the power, defined in W/mm, must be divided by the diffuser length.

Furthermore, the light source power was not considered as a parameter because it does not affect light propagation in tissues.

Results

a. Validation of the stereotactic registration

A phantom model was developed in order to evaluate the Leksell registration. A 15 cm 3 cube of polymethyl methacrylate (PMMA) was designed with three rods (see figure 4.6).

Thirteen radio-opaque marbles spaced 10 mm apart were inserted into these three rods. 10, 15, 20, 30, 40 and 50 mm)

c. Sensitivity study results

Sobol indices computed from the analytical model were illustrated according to the distance to the light source on figures 4.10,4.11 and 4.12. Figure 4.10 shows the first-order Sobol indices due to absorption coefficient variation. Figure 4.11 shows the first-order Sobol indices due to the reduced scattering coefficient variation. Figure 4.12 shows the second-order Sobol indices estimating the output sensitivity to absorption and reduced scattering coefficients. Table 4.4 summarizes the first-order and total-effect indices of each parameter (absorption and reduced-scattering coefficient). In this case, the total-effect index is the sum of the firstand second-order Sobol indices. These indices provide global sensitivity, which is the sensitivity of the parameter alone and interaction sensitivity with all other parameters. 4.4: First-order (𝑆 µ 𝑎 and (𝑆 µ 𝑠 ′ ) and total-effect Sobol indices (𝑆 𝑇𝑖 ).

Discussion

The use of a standard stereotactic procedure would facilitate the iPDT integration in a standard neurosurgery environment. Regarding the stereotactic registration algorithm implemented in the iPDT-Plan, its validation shows results close to the registration embedded in the Leksell GammaPlan® (Elekta, Stockholm, Sweden). The mean distance was 0.8 ± 0.3 mm, which is under the standard MRI resolution (i.e., one millimeter) commonly used in treatment planning.

As implemented in this study for µa = 0.02 mm -1 and µs' = 2 mm -1 , this analytical method (equation 4.9) shows an accurate estimation of light distribution, especially for a source length longer than 20 mm compared with the Monte-Carlo results. As detailed, the mean relative deviation of source length variations was approximately 1.23%, which is acceptable for a use as estimation of light propagation in a homogeneous medium, at a distance greater than the MFP' from the light source.

Regarding the sensitivity study, the reduced scattering coefficient is the most influential on analytical model variance close to the light source (approximately 70% to 80%).

Conversely, beyond a distance to the light source of 3 mm, the absorption coefficient becomes the most influential parameter on the analytical model. At a distance of 5 mm, approximately 98% of the analytical model variance is due to the effect of the reduced scattering coefficient.

First-order Sobol indices can be considered as total-effect indices between 0.75 mm and 10 mm from the light source, whereas the second-order Sobol indices are negligible. Between 10 mm and 20 mm from the light source, the impact of absorption coefficient on output variance remains constant, (approximately 20%) and the effect of reduced scattering decreases slowly (from 85% to 60%). Thus, the interaction between the absorption and the reduced scattering coefficient cannot be overlooked.

Since the analytical model implemented in the TPS does not consider the heterogeneity characteristic of GBM, the necessity to replace it by the Monte-Carlo model appears as highly recommended.

II. Forthcoming developments and improvements

Although this first attempt to design a TPS dedicated to the stereotactic iPDT of GBM enables keys steps to plan an entire iPDT procedure, several developments could be performed to improve the dosimetry part mainly, particularly to include the GPU Monte-Carlo model presented in Chapter III.

Automatic segmentation

A major advantage of the use of analytical models lies in its simplicity of implementation in a TPS. As seen previously, GBM is strongly heterogeneous, and the analytical models cannot consider this characteristic, which remains a major drawback.

As shown in Chapter III, the improvement of parallel GPU computing that drastically decreases the computing time [266,268,320] ensured strong enthusiasm for developing the Monte-Carlo method for PDT planning [102,266]. With the implementation of an automatic tissue segmentation methods [321,322] in the workflow of the iPDT-Plan (figure 4.1), the GPU Monte-Carlo model presented in chapter III could be used to compute the fluence rate according to optical coefficients for each voxel of the MRI volume used for planning.

Currently, the dosimetry model does not consider heterogeneity of target tissues. The entire space is considered as a homogeneous biological medium. A major improvement of the TPS will have to include an automatic segmentation of different parts of GBM, e.g., necrotic signal, enhancing and non-enhancing tumor and edema as proposed by the software BraTumIA (Brain Tumor Image Analysis) [START_REF] Velazquez | Fully automatic GBM segmentation in the TCGA-GBM dataset: Prognosis and correlation with VASARI features[END_REF]. This segmentation method, based on a decision forest classifier, modeled a map of the tumor cells distribution in tissue using four different MRI modalities: T1, T1-Gadolinium, T2 and FLAIR [START_REF]Traitements des gliomes de haut grade par thérapie photodynamique Les premières études de traitement de tumeurs cérébrales ont débuté en[END_REF]. From this map, distinction between four main biological tissue types is performed. Adding this fast automatic segmentation method in the actual workflow of the iPDT-Plan would enable to take advantage of Monte-Carlo possibilities, in particular to model light propagation in heterogeneous turbid medium. Thus, each voxel would represent a specific biological tissue with specific optical coefficient values (see table 4.5) at the segmentation output. Thus, the GBM heterogeneity would be considered. Presently, this software is not available in command-line interface and cannot be easily integrated in the iPDT-Plan workflow. Although, numerous algorithms can differentiate accurately different parts of brain structures, segmentation of brain tumors remains challenging, more particularly in the GBM case due to its strong heterogeneity [329].

Treatment optimization

The development of a computer-based optimization process occurring after the planning of cylindrical diffusers insertion would simplify the planning procedure. This also facilitates the standardization of PDT treatments, which enables to achieve better comparison between clinical studies.

Similarly to stereo-electroencephalography (SEEG) procedures [330,331]procedures which safety optimization have been recently proposed, iPDT planning for brain tumor therapy should include safety constraints during the optimization process. Indeed, SEEG planning tends to optimize the intracranial electrode implantation in order to maximize the recording volume, while constraining trajectories to safe paths. After the target points are defined, algorithm optimization computes possible trajectories and their associated risks.

Similar technologies may be used for iPDT. From fiber locations with associated risk along their path, light propagation in tissues might be computed to enable a global optimization process to illuminate a maximum of tissue volume while keeping trajectory safe.

Additional parameters integration

As previously developed in the second chapter (II, 2additional features for dosimetry planning), estimation of the distribution of light delivery is a major but not the only one component of PDT dosimetry. Photosensitizer and tissue oxygen concentrations should also be integrated in the workflow to obtain at least an explicit dosimetry of the PDT treatment planned. Thus, obtaining a complete dosimetry in clinic is a very challenging task which requires acquisition of several data before and during the surgical procedure with multiple probes insertion. Also, implicit features such as fluorescence and singlet oxygen assessment during the treatment could improve TPS dedicated to PDT. Presence of blood or cerebrospinal fluid can strongly affect the absorption of light during a PDT treatment. This alteration modifies the fluence deposited in tissue and, as a consequence, the PDT outcomes.

However, it remains difficult to prevent such perturbations on tissues optical properties and to integrate them into light distribution computing.

Finally, in analytical solutions, presence of non-diffusing part, numerical aperture and inhomogeneity of light distribution along the cylindrical diffusers are not considered. These features also influence the light emission. As Vesselov et al. showed in their study [332],

strong difference can appear between theoretical and experimental irradiation measurements from different manufacturers. Ideally, the efficiency and longitudinal, polar and azimuthal emission of each light source should be added into the fluence rate computation workflow to reach a more realistic model.

Chapter conclusion

A proof of concept of a TPS dedicated to the stereotactic iPDT is presented. This first attempt enables to plan the insertion of cylindrical diffusers and locate them in the Leksell frame, commonly used in neurosurgery. The stereotactic registration method used is described and its accuracy was evaluated by comparing its localization with the native commercial solution. Also, a simple dosimetry model based on analytical model, implemented in the software, enables to estimate light propagation in brain's patient. This model was compared to Monte-Carlo model to validate its accuracy. A sensitivity study was performed to evaluate the impact of absorption and scattering properties variations on the variance of the analytical model, as in GBM tissues.

Finally, several future developments are suggested, particularly to include a segmentation algorithm, which is mandatory for implementing the GPU Monte-Carlo model in the TPS workflow.

chapter presents a new lighting device dedicated to intraoperative GBM photodynamic treatment and for which the Monte-Carlo implementation presented in chapter III was deeply involved during the evaluation phase.

A short overview of previous clinical trials investigating intraoperative approach of PDT for GBM is first presented. Regarding the methodology expected to be seamless with the usual neurosurgery approach, a new lighting device was designed during this thesis and that includes specific expectations. Additionally, specific methodologies were developed to calibrate and assess the device in terms of mechanic and dosimetry. Indeed, the calibration lead to a transfer function associated with the device and that provides fast, robust and easy treatment duration prescription to induce a PDT response.

Furthermore, a comprehensive experimental design has been worked out prior to the clinical trial and results are reported in this chapter.

Finally, the setup of the clinical trial and first patient's results are introduced. Possible technological improvements and research perspectives are discussed.

I. Design and development of the prototype 1. Previous experiments

As previously described in the first chapter (II, 2 -Standard of care), when achievable, maximum tumor resection is expected according to the standard management of GBM [START_REF] Stupp | High-grade glioma: ESMO Clinical Practice Guidelines for diagnosis, treatment and follow-up[END_REF].

Fluorescence guided resection (FGR) leads to 2 months-increase of the progression free survival median (PFS). However, even if the surgery is complete (no more visible fluorescence or no more MRI signal highlighting the GBM), spread tumor cells usually remains within the healthy parenchyma. Because of these invading tumor cells, radiation therapy and concomitant and adjuvant chemotherapy fail to definitely cure the tumor and relapse is always occurring in the cavity margin in 85% of cases, [START_REF] Leroy | Fluorescence guided resection and glioblastoma in 2015: A review[END_REF]333,334]. A therapy to treat resection borders is highly expected to maximize the local control.

Illumination of the resection cavity remains challenging since the device has to suit the shape of the cavity. In that context, a deformable balloon is expected to meet this requirement.

As presented in chapter I (III, 3, a -Instrumentation), several devices have been designed in the past: inflatable balloon with different light sources [121-123, 132, 335], optical fiber coupled to a reflector in order to maximize light in the surgical cavity or LASER light source coupled to the surgical microscope and enabling a step by step illumination [124]. Moseley et al. described the only study using ALA as precursor [122,132]. For this study, both 5-ALA and Photofrin® were considered and repetitive PDT was evaluated. 5-ALA was used for the FGR procedure and PpIX was activated with the first illumination. The following illuminations (post-operative) were used to activate Photofrin®, which has a longer half-life than 5-ALA. The diffuser was created from a catheter balloon with a fiber and a cylindrical diffuser tip coupled to a Diomed PDT laser at 630 nm. Thus, none of these previous studies used 5-ALA as primary PS precursor to provoke a PDT effect on areas bordering the resection cavity.

Description of a new concept of intraoperative PDT procedure

To be more acceptable for both surgeon and patients, intraoperative PDT has to be fully embedded within the usual surgical workflow. Consequently, the concept of our new therapeutic approach was to deliver 5-ALA PDT intraoperatively early after FGR.

Accordingly, PDT added at the end of surgery achieved with FGR enables the treatment of the areas bordering the resection cavity while maintaining a maximal but safe resection (see figure 5.1). To make feasible intraoperative PDT, illumination of the surgical cavity must be carried out by a specific device diffusing a 635 nm LASER light homogenously to a sufficient irradiance.

Intraoperative PDT takes advantage of the presence of the PS which is already approved and administered to the patient for FGR purpose [START_REF] Leroy | Fluorescence guided resection and glioblastoma in 2015: A review[END_REF]. The proposed treatment strategy can fit into the current standard of care applied in first intention of care (see figure 5.

2) [START_REF] Stupp | High-grade glioma: ESMO Clinical Practice Guidelines for diagnosis, treatment and follow-up[END_REF]. Thus, intraoperative 5-ALA PDT is more ethically acceptable in comparison with interstitial PDT on resectable tumor. 2). This technic helps to better distinguish the healthy brain from remaining tumor tissue. In this case, the PS is used as a contrast agent for visual diagnosis.

The resection is ended when no remaining fluorescence is observed or to preserve functional areas. However, even if no fluorescence in the cavity is observed, remaining tumor cells around and deeper in the cavity are still present with PpIX uptake, but at a concentration limiting their visualization through the microscope. These remaining cells that are not macroscopically visible are the main target of the PDT treatment in that context. Thus, the PDT procedure is achieved early after the FGR. Here, PS is thus repurposed as a therapeutic agent to obtain a PDT effect. After the PDT procedure, as recommended by the standard of care, radio and chemotherapy are delivered to the patient to decrease recurrence risks.

To set-up an intraoperative PDT procedure, a lighting device has to be placed inside the surgery cavity and its inflation must be controlled by neurosurgeon until it reaches the cavity boundaries. The device must illuminate homogeneously the surrounding tissues with the help of an optical fiber. The treatment time enabling to achieve a therapeutic fluence value in surrounding brain tissues has to be easily and quickly computed.

For that purpose, a new lighting device was designed according to the comprehensive specifications drawn up in agreement with neurosurgeons from the Neurosurgery Department of the University Hospital of Lille. Additionally, a transfer function was built-up to enable real time dosimetry in the operating room according to the volume of the surgical cavity after gross tumor resection.

The medical device had to meet the following main requirements:

-Having a spherical or ovoid shape, -Owning a diameter up to 10 cm, -Enabling a conformability to the cavity, -Ensuring a biocompatibility with brain tissue, -Avoiding liquid injection in brain tissues, -Integrating a diffusing optical fiber at its center, -Being transparent to the light emitted by the LASER during the PDT treatment, -Allowing an easy setting up, -Enabling an acceptable illumination time.

Despite of the Herloon© balloon manufactured by B-Braun® Aesculap® (Tuttlingen, Germany) is neither dedicated to PDT nor neurosurgery, it meets several criteria of the specifications. Thus, to speed-up approval by ethical committees, this already CE marked and commercialized device was adapted for a PDT and served as basis for the design of the light applicator.

Medical device design

The device consists of two parts (see figure 5.3). The first one is the Herloon© balloon coupled to a trocar that is usually used during endoscopic intervention in which tissue needs to be separated in the extraperitoneal space. This device is built from a sterilizable trocar body and a single-use only transparent balloon. A sluice gate is added to the trocar, allowing fluid insertion to inflate the balloon. The inflation of the balloon ranges from a 4 cm diameter for 40 mL of diffusing solution to a 9 cm diameter for 500 mL of injected diffusing solution. A silicone valve terminates the trocar to ensure permeability and pressure inside the device. The second part, referred to as the optical fiber guide, was developed in the laboratory ONCO-THAI. This guide was built from an 8 mm diameter borosilicate glass tube closed at one end and inserted in a 10 mm stainless steel tube (see figure 5.4). This guide enables the positioning of the optical fiber at the balloon's center while ensuring complete permeability.

The optical fiber inserted in the guide is a 70 mm long cylindrical diffuser (RD-ML 70, Medlight, Ecublens, Switzerland) that approximately matches the balloon length in its deflated state. Finally, a screw thread fixed at the top of the guide allows tightening of the optical fiber inside the guide. The fiber guide is entirely compatible with the ethylene oxide sterilization process. To set up the whole lighting applicator, a single-use balloon is first inserted into the trocar. Then, the fiber guide hosting the 70 mm long cylindrical diffuser is inserted into the fiber guide (see figure 5.5) and is placed in the trocar balloon. After being assembled, the device is set on a frame and guided inside the surgical cavity by neurosurgeon. The balloon is filled with an emulsion composed of a 20% Intralipid® (Fresenius Kabi France, Sèvres,

France) diluted at a concentration of 0.1% until its wall reaches the brain-patient boundaries.

This emulsion exhibits weak absorption and strong diffusion at a laser light wavelength of 635 nm [122,132,335] and increases homogeneity of the light emitted from the source.

II. Calibration

The main purpose of the calibration was to model a transfer function so that the treatment duration can be predicted from an abacus according to the quantity of diffusing solution injected into the balloon. Indeed, the higher the quantity of liquid injected is, the larger the balloon diameter is and the longer the treatment duration has to be for a given therapeutic fluence. The treatment time, deduced from the abacus, is expected to lead to a therapeutic fluence of 25 J/cm² at 5 mm inside surrounding brain tissues.

The propagation of the light emitted from the balloon was quantified using several ex vivo experiments. To mimic clinical treatment conditions, calf brains were used due to their optical properties (absorption and reduced scattering coefficients of 0.19 cm -1 and 20.09 cm -1 , respectively [305]) close to those of human brains (absorption and reduced scattering coefficients of 0.20 cm -1 and 20 cm -1 , respectively) [117,336,337]. Additionally, these experiments were compared with in silico experiments from Monte-Carlo simulations using the algorithm implemented and detailed in chapter III.

Calibration of the measuring system

Before measuring fluence rate in calf brain using the same methodology and measurement chain than the Monte-Carlo validation experiments (III.2 Monte-Carlo validation by comparison with experimental measurements) calibration of the measuring system were achieved.

Indeed, as presented in chapter III (III, 2, a -Calibration factor), power measurements had to be converted to achieve the dosimetry of the device. Basically, a calibration factor CF 𝑏𝑟𝑎𝑖𝑛 was computed in order to convert power values measured in the calf brain P measured

[W] into fluence rate values ϕ brain [W.cm -2 ]:

ϕ brain = CF 𝑏𝑟𝑎𝑖𝑛 . P measured (5.1)
The same methodology as used in chapter III was applied using Marijnissen and Star papers [299,300]. Thus, four correction factors F 𝑛 , F 𝑝 , F 𝑏 , F 𝑖 and the calibration factor CF 𝑎𝑖𝑟 in air was estimated to obtain the calibration factor CF 𝑏𝑟𝑎𝑖𝑛 in calf brains.

ϕ brain = P measured . CF 𝑎𝑖𝑟 . F 𝑛 . F 𝑝 . F 𝑏 . F 𝑖 (5.2)
where F 𝑛 is a correction factor enabling to compensate the loss of light due to the interface sensor/surrounding calf brain using the refractive index of the calf brain around the isotropic probe (dimensionless). F 𝑝 is a perturbation correction [dimensionless]. F 𝑏 is a correction of the detection's photons' loss due to the surface occupied by the fiber on the detection sphere

[dimensionless]. F 𝑖 is a correction of the sensor's non-homogeneity response [dimensionless].

All factors were estimated using the experimental protocol already described in the third chapter (II, 2, a -Calibration factor).

Ex-vivo measurements

Prior to the experiments, fresh calf brains were prepared by dissecting out the cerebellum, the brain stem and the corpus callosum to improve the homogeneity of the studied tissues. The entire lighting applicator was fixed to a black box to avoid reflections and the balloon part was placed at the center of the box. Several pieces of calf brain were arranged in such a way that the balloon was entirely surrounded by tissues (see figure 5.6). To verify the validity of measurements, Monte-Carlo simulations were performed using balloon's dimension. The same protocol was used to compute treatment times according to the balloon's volume. Thus, the two transfer functions (ex-vivo and Monte-Carlo) were compared. Additionally, the difference between irradiance and fluence rate as previously introduced in chapter II is illustrated.

PpIX concentration impact on optical coefficients

It is known that a PS might induce a change in optical coefficients according to its concentration in tissues. Therefore, PpIX impact was evaluated using the equation described by .

The accumulation of PpIX does not affect the total reduced scattering coefficient. However,

III. Quality assessment and characterization

Several quality controls and characterization procedures were achieved in order to ensure the safety of the device. All these tests were used as technical reports for the setup of the clinical trial detailed later in this chapter.

Mechanical assessment

The behavior of the lighting device within a surgical cavity (pressure against the healthy parenchyma, conformity to the cavity) was studied by mean of a specific phantom with similar brain mechanical properties. Computerized Tomography (CT) images were acquired to evidence the conformity of the balloon when inflated in the surgical cavity.

a. Brain mechanical properties

The purpose was to evaluate the mechanical behavior of the device and its ability to fit into the surgery cavity. Two features can characterize mechanical properties in elastic domain: the Young's modulus and Poisson's ratio.

The Young's modulus describes the elastic behavior of a solid material. It represents the relationship between stress [N.m -2 or Pa] and strain [dimensionless] in continuum mechanics. A low Young's modulus (approximately kPa) describes a flexible material; contrariwise, a stiff material has a high Young's modulus (approximately GPa).

The Poisson's ratio describes the ratio of transverse strain to axial strain. Its value varies between 0 and 0.5. A high Poisson's ratio defines a material that can be deformed;

contrariwise, a material resistant to deformations or that breaks without showing any elongation owns a low Poisson's ratio.

Human brain has already been mechanically characterized, including Young's modulus and Poisson's ratio (see table 5.1).

Young's modulus (kPa) Poisson's ratio Miga et al., 2000 [340] 2.100 0.45 Miller et al., 2000 [341] 3.24 0.499 Ferrant et al., 2001 [342] 3.0 0.45 Four samples of different A:B ratio have been produced: (1:1.05), (1:1.08), (1:1.10) and (1:1.12) (see figure 5.9). Inside these four samples, three tensile specimens have been cut out. Twelve tensile specimens have been produced and characterized using extensometer.

This technic measures changes in the length of a material. Two points were drawn on each tensile specimen. By applying forces on each side of the tensile specimens, the two points Systemes, Vélizy-Villacoublay, France) for 3D printing in AM3300 polymer. Both, a head (see figure 5.11) and a waste mold of the brain negative were 3D printed. Using the suitable concentration mixing, the gel was poured into the waste mold. At the end, the resulting phantom gel owned exactly the same shape of the brain with the surgery cavity previously segmented. The phantom was imaged using CT-scan (see figure 5.12). The lighting applicator was placed inside the surgery cavity and filled by iodine solution in order to obtain a strong contrast between the balloon and the phantom's brain. Using this sequence, a conformation coefficient τ [%] was computed: τ = 𝐴 𝑐𝑜𝑛𝑡𝑎𝑐𝑡 𝐴 𝑐𝑎𝑣𝑖𝑡𝑦 (5.6) where 𝐴 𝑐𝑜𝑛𝑡𝑎𝑐𝑡 is the area where the balloon touches the brain's phantom and 𝐴 𝑐𝑎𝑣𝑖𝑡𝑦 is the area of the cavity.

Additional safety tests

In addition to the mechanical assessment of the device, several additional tests were performed to assess the safety of the device before being approved for a clinical trial.

a. Homogeneity of the light distribution

The homogeneity of the distribution along the balloon was evaluated. A flat sensor (818-SL, Newport, Irvine, CA, USA) was placed against the balloon at five different locations (see figure 5.13). The irradiance was measured according to five different quantities of diffusing solution injected into the device (100, 200, 300, 400 and 500mL). For each volume, the irradiance was normalized.

c. Thermal variations

The temperature elevation was also investigated to evaluate the temperature rising risk of surrounding brain tissues (see figure 5.14). The balloon was filled with 500 mL of diffusing solution and then immersed in a hot-water bath adjusted to 37°C. The LASER has emitted a continuous power of 2 W for two hours. The temperature was measured with a thermographic camera. IV. Results

Dosimetry a. Calibration factor in calf brains

The calibration factor CF 𝑏𝑟𝑎𝑖𝑛 was estimated using the following values of correction factors F 𝑛 , F 𝑏 , F 𝑝 and F 𝑖 (see table 5.2). The correction factors F 𝑛 , F 𝑏 and F 𝑝 were determined from the papers of Marijnissen et al. [299,300]. The non-homogeneity response factor F 𝑖 was determined from the standard deviation of measurements where the isotropic probe is irradiated according several irradiation angle θ (figure 3 This calibration factor value is only valid for the isotropic probe used during the power measurements and calf brains. Relative error between fluence rate values against the balloon's wall obtained from power measurements and Monte-Carlo method were computed (see table 5.4). By fixing a fluence value of 200 J/cm² against the balloon's wall, a transfer function has been generated from the Monte-Carlo results (see figure 5.19). The relative error between the two transfer functions equals to 9.998%.

Volume of the

d. PpIX impact on optical coefficients

In the study of Valdes et al. [352], a PpIX concentration of 4.523 µg/mL (approximately 8.039 mol.L -1 ) was observed in GBM tissue with a high level of fluorescence and PpIX accumulation. Using this concentration and at 635 nm, the PpIX absorption coefficient was estimated to be 0.0187 cm -1 in a way that the total absorption coefficient was only increased of 8.5% (µ 𝑎,𝑡𝑜𝑡𝑎𝑙 = 0.2187 𝑐𝑚 -1 ). PpIX impact on optical properties is negligible at 635 nm.

Mechanical tests results

a. Evaluation of the phantom Young's modulus Figure 5.21 shows the curves obtained from the eleven tensile specimens characterization using extensometer: 3 tensile specimens of (1:1.12), (1:1.10), (1:1.08) and two tensile specimens of (1:1.05) (one (1:1.05) tensile specimen has teared during the test).

The same stretching had been applied on each tensile specimen. Nominal stress [START_REF] Kawczyk-Krupka | ALA-induced photodynamic effect on vitality, apoptosis, and secretion of vascular endothelial growth factor (VEGF) by colon cancer cells in normoxic environment in vitro[END_REF] is plotted according to the stretch [mm]. Modeling using Yeoh's model is superimpose on each acquisition. Young's modulus has been computed using Yeoh's model for each mixing (see table 5.5).

Mixing proportion of B (%)

(1:1.12) (1:1.10) (1:1.08) (1:1.05)

Young's modulus of tensile specimen 1 (kPa)

c. Thermal variations

Regarding the temperature measurements, no temperature elevation has been observed

during the experiments. Thus, the temperature elevation due to the device during treatment can be considered acceptable.

d. Robustness

Regarding to the robustness evaluation, a maximal volume of 1600 mL of diffusing solution has been injected into the device until balloon breakage.

V. Clinical evaluation 1. Setting up of the clinical trial

Once the lighting applicator was characterized and validated, a phase one clinical trial has been investigated to evaluate the innocuousness of the procedure. This key milestone in the development of a medical device required several studies and particular methodology.

Thus, the pilot study has been setting up to assess feasibility and safety of intraoperative 5-ALA PDT procedure for the treatment of newly diagnosed GBM. This phase one clinical trial is monocentric, prospective is leaded by Pr. Reyns as principal investigator and Dr. Vermandel as scientific coordinator.

The primary endpoint is to assess the feasibility of intraoperative PDT early after the surgical resection of GBM without unacceptable and unexpected toxicities.

The secondary endpoints are the evaluation of:

-progression-free survival (PFS),

-overall survival (OS),

-treatment response on MRI, every three months until tumor relapse, -toxicity every three months until tumor relapse and follow-up of serious unexpected event potentially linked to patient's death, -quality of life, every three months until tumor relapse.

Enrollment of 10 patients is expected over a period of 12 months. Approximately 150 new cases of GBM are diagnosed each year in the Lille Hospital. Complete resection is usually feasible for more than a half of these cases. Evaluation criteria have been established in order to evaluate objectively the results obtained. The main criterion that enables the validation of primary endpoint, is the evaluation of the patients' proportion receiving full PDT treatment (expected to be at least 70%), according to the following requirements:

-achievement of the full PDT treatment with an acceptable tolerance due to the -Assessment of quality of life using a questionnaire in pre-treatment and posttreatment at 3, 6 and 12 months after treatment during neurosurgery consult, -Treatment response on MRI at 3, 6 and 12 months after treatment.

A user manual has been written describing the setting procedure from the device assembly until the end to the treatment. A technical report has been edited describing all different parts of the device; characterization and safety tests were realized to insure its harmlessness. Materials and their biocompatibility have been defined. Ex-vivo experiments and sterilization process were also presented. Risks analysis was achieved in accordance with the European directive 93/42 EEC on medical devices and according to the European standard EN ISO 14971: 2007. Other reports had to be written to complete the application to authorities, such as investigator brochure, case report form or record sheets to organize the quality control system.

All these documents were submitted to the two French competent authorities, e.g., the ANSM ("Agence Nationale de Sécurité du Médicament et des Produits de Santé") and the ethical committee CPP ("Comité de Protection des Personnes") for technical (benefits and risks) and ethical aspects approvals respectively.

First patients

On September 22nd and December 26th of 2016, the CPP and ANSM approved the setting of the clinical trial. The first patient has been enrolled on 5th May of 2017. In the meantime, a patient test using 3D printed patient's head harboring a surgery cavity had been programmed to achieve a complete procedure with the neurosurgery team and rectify issues that could have appeared (see figure 5.24). Several features of the PDT procedures have been measured (see table 5.6). This study demonstrated that a value of 0.2 cm -1 appears to be an absorption coefficient close to mean human brain optical properties, without considering gliomas cells infiltration [117,305,336]. Based on this value, modeling of the function transfer with Monte-Carlo is close to the one based on ex vivo experiments. This similarity reinforces the validity of the transfer function applied to the intraoperative PDT. However, ex vivo experiments could not be performed after 5-ALA administration, and modelling PpIX metabolism was not achievable. Thus, the impact of the absence of PpIX in the biological sample was evaluated prior to experiments. Indeed, differences in the absorption coefficients with and without PpIX might lead to different fluence rates. Nevertheless, the PpIX absorption coefficient in this condition and at 635 nm is only responsible for an increase of 8.5% in the total absorption coefficient, and the impact of PpIX was thus considered negligible for ex-vivo experiments. Moreover, the tissue surrounding the balloon would have a lower PpIX uptake than resected tumor tissues and the corresponding PpIX impact should still be more modest, confirming the assumption that the ex vivo model used in this study suited well the characterization of the device.

Patient

Regarding the precision of converting power measurements into fluence rate values, the calibration factor requires several constants including the refraction index of the medium surrounding the probe. This data strongly increases the calibration factor estimated in the air.

The impact of the standard deviation of the calibration factor in the biological medium induces a standard deviation of 62 seconds for the illumination time (approximately 11% of the treatment duration for a mean cavity volume considered in this study). This value remains low compared with the total illumination time needed to treat one patient. A slight influence of calf brain oxidation on the power measurements, which cannot be quantified, was observed. However, to minimize this variable, fresh specimens preserved by refrigeration

Clinical feedback

Although several studies have already reported intraoperative PDT, some experimental conditions have changed since the last trial, and more particularly the PS administered to the patient. In the last clinical trial in 2007 [122,132], the drug administered was Photofrin ® (Pinnacle Biologics™, Chicago, IL, United States). This first-generation PS had some disadvantages, including low absorbance at 630 nm and poor selectivity for tumor cells, leading to the necrosis of healthy tissues [354]. With the new generation of PSs, the selectivity and reactivity have been improved. In the surgical procedure investigated, the precursor 5-ALA (Gliolan-Medac, Germany) is administered to the patient, which induces a strong PpIX uptake in tumor cells [39,355]. With the application of the device soon after FGR, improvement in the local control of direct PDT effects and an immune response with indirect PDT effects are expected. Including the 5-ALA PDT inside the SOC makes this trial a world premiere.

Additionally, investigating the impact of PDT on the patient's immune system is highly expected. Indeed, exosomes [356], which may be potential biomarkers of the PDT effect, can be massively released spontaneously by tumors after PDT treatment. Observing its evolution in the course of the management of the patient, the efficacy of PDT treatment could be related to the amount of exosomes. Such an analysis would enable to carry out immunemonitoring of patients by controlling the various immune populations (e.g., T lymphocytes, B lymphocytes, Regulatory T cell). Medlight, Ecublens, Switzerland) during the illumination would allow to record the light inside brain tissues [191,[194][195][196][197][198]. Spectroscopy technic can detect the fluorescence peak and monitor its intensity during the procedure. Also, by modifying the depth insertion of the probe, optical properties of surrounding tissues can be estimated. Regarding the singlet oxygen monitoring, a direct measurement of the very weak luminescence peak in the near infrared can only be performed using the right equipment.

VII. Forthcoming improvements

b. "Smart balloon": photo-acoustic imaging

Recent studies show an interest in the photo-acoustic technology to imaging changes in biological tissues [357,358]. This imaging modality is based on the emission of ultrasonic waves by biological tissues irradiated. When laser beam is absorbed by a tissue, this energy absorption produces a local heat increase. A fast thermo-elastic expansion is created, leading to a displacement of particles that causes a pressure wave. By recording these waves with transducers, it is possible to determine the localization of the wave's origin and the type of absorbing material. With the miniaturization of transducer technology, capacitive micromachined ultrasonic transducers (CMUT), which can reach a size of 10 to 100 µm, could be placed around the surface of the balloon to perform a real-time monitoring of the PDT treatment. One of the main features monitored is the blood oxygen saturation.

Photoacoustic tomography is used to record vascular damage due to PDT treatment in a high resolution (approximately 100 µm or less). However, to obtain such data, additional laser sources are necessary to excite molecules targeted: 750 and 850 nm for the study of Mallidi et al. [359] or 560 and 576 nm for the study of Shao et al. [358]. Although this technology can provide depth imaging of PDT impact during the illumination, it still remains far from being used in clinic.

Lighting fabrics

Another perspective of improvement lies in the use of light emitting fabric for intraoperative PDT. This technology already evaluated in dermatology clinical trial [360,361] could be adapted to intracavitarly PDT treatment. The term of "light blanket" has already been used in intracavitarly context [362][363][364]. However, the technology was different from the actual light emitting fabric. The light blanket was composed of a series of parallel cylindrical diffusers, fixed with a polyethylene blanket and two other layers (a diffusing solution compartment and a reflector). The light emitting fabric uses the Fresnel conditions to emit light. If an optical fiber is bended beyond a critical angle 𝜃 𝑐 (Brewster's angle which depend on the refractive indices of optical fiber materials), the light leaves the optical fiber:

𝜃 𝑐 = arcsin ( 𝑛 2 𝑛 1 ) (5.8)
where 𝑛 1 and 𝑛 2 are respectively the refractive indices of the fiber core and surrounding cladding.

Using different weaving technics of optical fibers optimized for the light leak, called side emitting optical fiber, lighting fabric of different surface areas can be obtained. The main advantage of this technology is its easy conformation to all different topologies of treatment surface, including intracavitarly treatment. Nevertheless, current light emitting fabrics used in dermatology reach an irradiance of approximately 18 mW/cm². To perform a PDT treatment with such a system in order to replace the actual proposed balloon device, optimization had to be performed to increase the irradiance delivered. Also, in clinical context, these fabrics would have to be sterilized. Strict restrictions are applied with class III medical device (medical device in direct contact with central nervous system). Although this alternative seems to be interesting, several challenging studies had to be achieved to evaluate it in a clinical trial dedicated to intraoperative PDT GBM treatment.

Chapter conclusion

Unfortunately, the PDT effect on GBM is still not evidenced because of the poor reproducibility of the procedure from one center to another, the multi-modality treatment when considering only the recurrent GBM population and the low inclusion numbers in the past trials.

In this chapter, a new device and the associated procedure is proposed in order to facilitate the incorporation of intraoperative PDT into the standard surgical workflow. The innovation is better assimilated when its users are involved in its design for they adapt it to their expectations. Thus, the design of this new device aimed to simplify the PDT procedure for the neurosurgeon with a seamless methodology suitable for the standard surgical workflow which also limited the dosimetry issue during surgery. In such a manner, and thanks to a highly reproducible procedure, the effects of PDT should be better highlighted on a larger and more homogeneous population of only newly diagnosed GBM patients.

Although the device presented here could be improved, in particular with its adaptability to the cavity, this lighting delivery system is a first step towards future larger clinical trials and a future integration in the GBM standard of care.

General conclusion

The main purpose of this thesis work was to provide technologies to enable the GBM treatment with PDT. In order to include the totality of patients harboring GBM (resectable and unresectable cases), two different guidelines have been elaborated.

For unresectable tumors, the stereotactic iPDT is proposed. Because such a surgical procedure had to be planned, a TPS dedicated to this technic must have been developed.

Currently, no TPS dedicated to iPDT for GBM treatment is available. In Chapter IV, a proof of concept of such a software is presented. Although several improvements are needed, it enables both the localization of cylindrical diffusers in the stereotactic Leksell frame usually used in neurosurgery and the estimation of the light propagation emitted from these light sources without considering the heterogeneity of GBM tissues. In order to improve this estimation, Monte-Carlo method, often used as reference to compute the light distribution, should be include into the workflow. However, because the main drawback of this model lies in the large time required to perform an entire simulation, its use in clinical routine remains an issue.

The improvement of parallel computing using GPU that drastically decreases the computing time [266,268,320] Regarding the performances, the reported acceleration (approximately 660 times) is in accordance with the literature. The time required to compute the light propagated from a cylindrical diffuser is now sufficiently fast to be integrated in clinical routine.

For resectable GBM cases, an intraoperative PDT modality treatment is proposed early after FGR. In Chapter V, a new device dedicated and its inclusion into the standard of care are described. The main asset of the proposed surgical procedure is to provide a straightforward use of the device and a surgical procedure suitable for the standard surgical workflow. The This device can fit into large surgical bed (>150cm 3 ). However, as the mechanical assessment and clinical evaluation show, when the device is placed in a small cavity, the balloon is not completely inserted. The impact on the deposit light during the treatment is negligible and a PDT effect is still expected at margins resection. To circumvent this issue, other devices such as light emitting fabrics should be assessed.

As previously announced in the introduction, the essential goal of research in medicine is patient's treatment. Accordingly, the feasibility of the delivering 5-ALA PDT for the treatment of newly diagnosed GBM is currently under evaluation in the world premiere INDYGO clinical trial (ClinicalTrials.gov identifier: NCT03048240). Ten patients will be recruited in this phase one clinical trial.

To conclude, this thesis moves the PDT treatment of GBM to a new step. Although previous results show inspiring outcomes, including some long-sustaining response patients, the effects of PDT should be better highlighted on a larger and more homogeneous population of newly diagnosed GBM patients only.

A strong and multidisciplinary effort should be conducted to propose a combined therapeutic solution, particularly with PDT and immunotherapies, in order to find a cure for brain cancer. Plusieurs sources lumineuses sont développées (fibres diffusantes, ballonnet diffusant, tissus lumineux) pour adapter l'illumination à la cible thérapeutique (voir figure 1.7).

List of figures

Le glioblastome

Les gliomes sont les tumeurs du cerveau les plus fréquentes chez l'adulte. Elles sont classées en fonction de quatre grades par l'Organisation Mondiale de la Santé [4]. Le glioblastome (GBM) ou gliome de grade 4 est la tumeur la plus fréquente et la plus agressive interstitielle (iPDT) peut être envisagée. Cette modalité de traitement consiste en l'insertion de fibres optiques (diffuseurs frontaux ou cylindriques) au sein même de la cible thérapeutique afin de délivrer la lumière directement in situ. Dans la plupart des cas de traitement de GBM par iPDT, les fibres optiques sont insérées dans un référentiel stéréotaxique qui permet de localiser leur insertion dans le cerveau du/de la patient.e. Les principaux essais cliniques de traitement de GBM par iPDT sont listés dans le tableau 1.1.

Dans le cas où la chirurgie peut avoir lieu, la PDT peropératoire (perPDT) peut être administrée en fin de résection. Un dispositif lumineux (ballon, fibre cylindrique, laser couplé au microscope chirurgical) est inséré dans la cavité afin de générer un effet cytotoxique dans les berges opératoires. Les principaux essais cliniques de traitement de GBM par perPDT sont listés dans le tableau 1.2.

Dosimétrie prédictive de la PDT

La dosimétrie prédictive utilisée pour planifier la délivrance d'un traitement PDT comprend la modélisation des trois composants principaux de la PDT (i.e., la propagation de la lumière, la distribution du PS et la concentration de l'oxygène dans les tissus).

Ainsi, la modélisation du transport de la lumière dans la matière est une composante des cellules de gliomes [14,[30][31][32].

La puissance [W] est l'unique grandeur pouvant être mesurée pour quantifier la quantité de photons en un point de l'espace. Cependant, pour pouvoir calculer une dose (ou fluence [J/cm²]), l'estimation du débit de fluence [W/cm²] est nécessaire. Un facteur de calibration permettant de convertir les mesures de puissance en débit de fluence a été calculé à partir des études de Star et Marijnissen [33,34], recommandées par l'American Association of Physicists in Medicine (AAPM) [35,36]. Ce facteur de calibration est le résultat de la multiplication de cinq facteurs permettant de corriger l'influence du capteur sur la mesure (inhomogénéité, dimensions du capteur, matériaux utilisés).

Une fois ce facteur correctif estimé, cinq fibres cylindriques ont été placées au centre sensibilité a été évaluée à l'aide des indices de Sobol [37,38]. Les résultats de cette étude montrent que l'absorption pondère plus la diffusion à une distance supérieure de 3 mm de la source lumineuse, et inversement pour la diffusion. Cette sensibilité de la variance de la méthode analytique souligne l'importance de remplacer le modèle analytique par le modèle Monte-Carlo. présentes dans les berges de la cavité. Ce résidu de cellules est à l'origine de plus de 85% de récidive tumorale [15,39,40]. Un contrôle local pour diminuer ce risque de récidive est donc primordial. La PDT intervient donc en fin de FGR, profitant ainsi de l'administration du PS au/à la patient.e pour le fluoroguidage. Le PS est ainsi repositionné an tant qu'agent thérapeutique et non de contraste. Pour réaliser cette illumination, un dispositif dédié à la perPDT a été développé.

Ce dernier repose sur les études antérieures réalisées pour développer l'illumination peropératoire de la cavité opératoire [41][42][43][44][45][46]. Le dispositif présenté dans cette étude est composé de deux parties. La première, commercialisée par B-Braun® Aesculap® (Tuttlingen,

  therapy ................................................................................................... 23 1. LASER: History ................................................................................................... 23 2. Photosensitizer ..................................................................................................... 28 3. Photodynamic therapy .......................................................................................... 33 II. Glioblastoma ................................................................................................................... 38 1. Description and latest classification ..................................................................... 39 III. PDT treatment of HGG ............................................................................................. 43 1. History .................................................................................................................. 43 2. Interstitial PDT ..................................................................................................... 44 3. Intraoperative PDT ............................................................................................... 46 Chapter conclusion ................................................................................................................. 50

Figure 1 . 1 :

 11 Figure 1.1: Illustration of the stimulated emission of radiation.

Figure 1 . 2 :

 12 Figure 1.2: Jablonski diagram of triple state inversion population.

Figure 1 . 3 :

 13 Figure 1.3: Laser evolution from 80's to current device.

Figure 1 . 4 :

 14 Figure 1.4: Jablonski diagram of photophysical processes.

  nor heme is accumulated in a normal cell. In the case of tumor cells, individual steps of the heme synthesis process are altered and might cause accumulation of various intermediate products leading to a PpIX accumulation in tumor cells only. Then, an exogenous exposure to 5-ALA leads to a strong escalation of PpIX concentration in tumors cells up to 6-8 h before saturation.

Figure 1 .

 1 Figure 1.5: a) 5-ALA molecular structure, precursor of the photosensitize b) Protoporphyrin IX molecular structure.

Figure 1 . 6 :

 16 Figure 1.6: Illustration of the light depth penetration in biological tissue: a) dependency of the wavelength and b) dependency of the size of the light source (taken from[36]).

Figure 1 . 7 :

 17 Figure 1.7: Illustration of lighting modalities available for PDT treatments (illustrations taken from [41]).

Figure 1 . 10 :

 110 Figure 1.10: Decision algorithm of the HGG treatment plan.

[

  101], urology [102-104]. Some more advanced or commercial applications are available such as iDOSE developed by Spectracure [105]. It provides an optimization of the fiber position into the body's patient, generated from a 3D tissues model using ultrasound imaging. A second example is the TPS provided with the technical solution named Tookad developed to plan an interstitial PDT treatment of prostate cancer [106].

Figure 1 . 11 :

 111 Figure 1.11 : The two main clinical PDT modalities for neurosurgery a) Interstitial PDT: optical fibers are inserted through the skull into the tumor core that contains heterogeneous tissues (necrosis, tumor, and edema). b) Intraoperative PDT: a balloon device is inserted into brain patient to treat borders of the resection cavity.

Figure 2 . 1 :

 21 Figure 2.1: Light absorption process by a chromophore of a geometrical cross-section A [m 2 ] and an effective cross-section σa [m 2 ] linked by the absorption efficiency Qa [dimensionless].

Figure 2 . 2 :

 22 Figure 2.2: Harmonic oscillator modeling the scattering process.

Figure 2 . 3 :

 23 Figure 2.3: Relation between wavelength and associated scattering theory.

Figure 2 . 4 :

 24 Figure 2.4: Scattering process by a chromophore.

Figure 2 . 5 :

 25 Figure 2.5: Henyey-Greenstein functions for different anisotropic values: plot a) and polar plot b) of the scattering phase function p(θ) [sr -1 ] according to the angle θ (The forward direction along the original photon trajectory is 0° and the backward direction is 180°).

Figure 2 . 6 :Figure 2 . 7 :

 2627 Figure 2.6: Mean free path versus reduced mean free path. The term "reduced" enables to extrapolate several progressions of one photon into one mean displacement.

Figure 2 . 8 :

 28 Figure 2.8: Radiance illustrated.

Figure 2 . 9 :

 29 Figure 2.9: Irradiance computing illustrated. A light beam is emitted on a surface; the irradiance on the area containing the entire beam is computed.

Figure 2 . 10 :Figure 2 . 11 :

 210211 Figure 2.10: Irradiance versus fluence rate: in A), the fluence rate Φ is three times greater.than the irradiance E. In B), the irradiance E is defined as the power striking the surface.Conversely, the fluence rate Φ considers all photons distribution crossing the cross-sectional area. Due to the high scattering property of biological media, the fluence rate value is higher than the irradiance value.

Figure 2 . 12 :

 212 Figure 2.12: Analytical models estimate the fluence rate radiated from different types of light sources : A) light point source (equation 2.25) B) cylindrical light source considered as sum

  The three main element methods allowing to solve the diffusion equation are the finite difference, finite element and boundary element methods. All these methods mostly differ from the initial geometry of the mesh (see figure2.13). The finite difference method discretizes the entire space into a Cartesian mesh. The finite element discretizes the region of interest in 2D with triangular elements or in 3D with tetrahedral elements. The boundary element discretizes boundaries of the different homogeneous parts of the region of interest in several nodes.

Figure 2 .

 2 Figure 2.13: A) Cartesian mesh used in finite difference method B) 3D mesh with tetrahedral elements used in finite element method C) 3D mesh with several nodes used in boundary element method.

Figure 2 . 14 :

 214 Figure 2.14: Illustration from Comsol MultiPhysics simulation using finite elements method to resolve different PDE. Here, the Helmholtz equation is solved: A) two cylindrical diffusers are defined as light source and B) an adaptive mesh is built around these cylinders in order to obtain an accurate results near boundaries (C and D).

  photon radiation therapies. Wilson and Adam performed the first Monte-Carlo simulations of light (non-ionizing photons) transportation in a biological medium in 1983 [166]. Since then, Wang, Jacques and Prahl [167, 168] considerably improved Monte-Carlo performances and accuracy leading to several Monte-Carlo codes including the well-known MCML program (Monte-Carlo Multi-Layered) [169].

Figure 2 . 15 :

 215 Figure 2.15: Illustration of the Monte-Carlo method. The area of the lake is determined from the number of balls fired randomly from a cannon fallen or not in the lake. Higher the number of balls fired is, better the estimation of the lake area is.

2 ]

 2 [S0], [S1] and [T] are respectively the singlet and triplet PS concentration, and are respectively the ground triplet and excited singlet state oxygen concentration, P and [A] are the oxygen supply term and the concentration of [ 𝑂 1 2 ] acceptor excluding PS molecule and 𝛿 is a low PS concentration correction term. The ki parameters define different rates (photon absorption, decay excited to ground state or bimolecular rate of reaction) and are described in the Wang et al. study [180]. These six differential equations (equations 2.33) are then combined with the light diffusion equation and solved in order to obtain the concentration of the excited singlet state oxygen[ 𝑂 1 2 ], which induces the cytotoxic effects of PDT treatments. In the study of Liang et al., this method was implemented in the software Comsol Mutliphysics® dedicated to numerical analysis using finite element method [182]. This approach was evaluated on C3H mice treated with two PS: Photofrin and BPD Verteporfin. The authors also associated the concentration of excited singlet state oxygen to the size of necrotic areas predicted and measured. More recently, Qiu et al. correlated three parameters: the fluence rate, the fluence and the concentration of excited singlet state oxygen [183]. The authors argued that the parameter [ 𝑂 1 2 ] was the most accurate metric and could serve as a better predictive dosimetric quantity. Additionally, Lopez et al. described a new dosimetry metric: the tumor reactive single oxygen as a new dosimetry metric that represents the amount of singlet oxygen per tumor volume that reacts with the molecules in the tumor [184]. The authors proposed in their study a model using the PPIX as PS in 5-ALA PDT treatment. This metric integrates fluence rate, optical properties, the oxygen consumption and the PS distribution.

  : from basics to PDT implementation I. Monte-Carlo model .......................................................................................................

  ........................................................................................................... 103 3. Monte-Carlo validation by comparison with experimental measurements: results108 4. Performances ...................................................................................................... 111 Discussion and chapter conclusion ..................................................................................... 116

Figure 3 . 1 :

 31 Figure 3.1: A photon history; from launch to termination.

  photons simulated, Nstep is the number of steps required for each photon and 𝑡ℎ𝑟𝑒𝑠ℎ𝑜𝑙𝑑 is typically 10 -4 . Those factors are computed using previous Monte-Carlo simulations and are valid only for a given computing configuration.

Figure 3 . 2 :

 32 Figure 3.2: Flow-diagram of a standard steady-state Monte-Carlo simulation.

  had access to the common disk during one hour per day. However, with the advent of multiple Central Processing Unit (CPU) or multiprocessing, new possibilities to speed up Monte-Carlo simulations have emerged. Colasanti et al. developed a code addressed to multiple CPU contained into one unit computer. The high performances exhibited by this new type of processor (cost about 2500$) considerably speed-up the resolution of light transport [256]. Reported acceleration reduced the computing time by a factor equal to the number of CPU used. It ranged from a factor 100 for a number of 128 CPU to a factor 1000 for a number of 1024 processors, the maximum capacity available in 2000. A recent coprocessor (Intel Xeon Phi) was used to simulate light propagation in turbid media [257].ii. Internet-based parallelismThe Internet revolution has also strongly impacted performances of high computational algorithm. This technology recently used in treatment planning systems mostly dedicated to radiotherapy tends towards its application into the simulation of visible photons simulation[258][259][260][261]. Thus, using the mapReduce3 environment provided by Google to facilitate development of parallel algorithm including cloud computing, Pratx et al. created a cloud version of the MC321 Monte-Carlo code by Steven Jacques [262, 263]. This monte-Carlo code simulated the photons propagation in steady state in a homogeneous and infinite medium and results could be stored in spherical, cylindrical and planar coordinates. Performances of the fastest approach described in the paper reached 2 minutes for simulating 100.10 6 photons, compared to 28 minutes for a single thread. A speed-up of 1258 was reported in comparison with a value of 1080 with the fastest Monte-Carlo method available in 2011 on a single computer. A peer-to-peer Monte-Carlo code was also developed by Doronin et al. in order to simulate photon transport in heterogeneous turbid media [264]. c. GPU acceleration Recently, one major advancement in parallel computation was realized using the Graphics Processing Unit (GPU) to speed-up considerably Monte-Carlo simulations at a limited cost. The first article introducing this new technology was the study of Erik et al. in 2008 showing that a low-cost GPU (110$) could provide extreme acceleration [265]. Authors transposed a code used to simulate photons migration in a semi-infinite, homogenously scattering and non-absorbing medium in a Compute Unified Device Architecture (CUDA) architecture [224]. This platform, ownership of the trademark NVIDIA (NVIDIA Corporation, Santa Clara, CA, USA), was released in 2007 and enabled to develop a C, C++ or Fortran code in order to take advantage of the computing parallel possibilities offered by a GPU. To compare the performances of their code, authors simulated 500,000 photons in a White Monte-Carlo model (i.e., a Monte-Carlo simulation of a homogeneous medium with zero absorption) on both CPU and GPU. The speed-up was about 1000 times faster with the GPU: 8513 seconds for the CPU, 7.9 seconds for the GPU. Their work was enhanced two years later by proposing a GPU version of the MCML code [169] and optimizing fast access memories. They reported a GPU acceleration of about 600 times faster than CPU version in heterogeneous media. The most popular GPU version of time-resolved photon propagation in a 3D turbid media was proposed by fang et al. in 2009 [266]. They claimed an improvement of 300 times faster than conventional CPU simulation. Other simulations were then developed using GPU parallelization. Martinsen et al. implemented a photons transport simulation through turbid media and reported an acceleration of 70 times faster than CPU only [267]. Several other studies also proposed to simulate photons through turbid media with several discrepancies [268-270]. Carbone et al. introduced the GPU computation in order to accelerate the solving of the forward problem to retrieve optical coefficients of inhomogeneities introduced in homogeneous medium [271]. Leung et al. simulated ultrasound-modulated light in turbid media with an improvement of 125 times of the standard computing time [272]. Powell et al. took the advantage of GPU acceleration to simulate acousto-optic effect in heterogeneous turbid media [273]. Yang et al. generated a diffuse reflectance value in 0.08ms in order to interpret quickly fiber-based and camera-based measurements of remitted light from biological tissues [274]. The more recent article, provided by Li et al. compared the performances of 3 different GPU in order to study the behavior of polarized photons as they propagate through complex biological tissues [275]. All these papers demonstrate a strong interest in this technology that provides the better costefficient choice to improve Monte-Carlo simulations performances.

)

  From a sequential point of view, this sum requires 𝑚 × 𝑛 consecutive sums. If a sum takes tsum time, this operation necessitates 𝑡 𝑠𝑢𝑚 × 𝑚 × 𝑛 without taking into account access time and reading/writing time on the hard disk drive that contains values of matrices. However, from a parallelized point of view, because all 𝑚 × 𝑛 sums are performed simultaneously, the entire computation takes only 𝑡 𝑠𝑢𝑚 time. This example easily shows the gain provided by such implementation (see figure3.3).

Figure 3 . 3 :

 33 Figure 3.3: Sequential versus parallelized computing: time required to perform the sum of two 𝑚 × 𝑛 matrices is decreased of a 𝑚 × 𝑛 factor.

Figure 3 . 4 :

 34 Figure 3.4: Each kernel hosts the parallel part of the code; several kernels are gathered into a block and all blocks constitute a grid.

Figure 3 . 5 :

 35 Figure 3.5: The six different memories constituting the device memory and their read and write access.

  computing platform Numerous parallel computing platforms have been developed to communicate with the device. Among them, two are commonly used to develop applications on GPGPU: Compute Unified Device Architecture (CUDA) and Open Computing Language (OpenCL).In 2007, the American company NVIDIA, mostly known for its GPU production, released the CUDA platform. This parallel computing platform and application programming interface (API) was developed to facilitate the communication with GPU and introduce the use of kernels as processing units. CUDA was elaborated with programming languages such as C, C++, and Fortran. This platform can only work with the use of NVIDIA GPU. In return, NVIDIA provides a well-integrated library allowing an easier control of GPU programming.In 2008, Apple Inc. developed a proposal of framework dedicated to GPGPU. In 2009, the Khronos Group, a non-profit technology consortium, released the initial version of this framework. OpenCL was elaborated with C and C++ programming languages. The main difference with its CUDA rival lies in the opening to other GPU brands such as AMD. The use of OpenCL does not require a NVIDIA GPU. However, codes developed on this platform will not run optimally on all of them without strong optimization work.b. ImplementationOur Monte-Carlo method has been designed toward its implementation in a future treatment planning system (TPS) dedicated to photodynamic therapy for glioblastoma (GBM) treatment as introduced previously in chapter I. Monte-Carlo simulations rely on segmented MRI examinations for which each voxel is classified according to the type of biological medium it belongs to (white matter, grey matter, tumor, edema, necrosis). This data represents a 3D matrix in which light sources are modelled according their location and shape, basically one or several cylindrical diffusers.Structures ofFang et al. [266] and the mcxyz CPU code by 279]. Several input features are required to execute the simulation. 12 parameters are written in a first text file:-Number of voxels on each axis (Nx, Ny, Nz), -Voxels size on each axis (dx, dy, dz), -Coordinates of a first point, localized on the x, y and z-axis, corresponding to the beginning of the cylindrical diffuser (xin, yin, zin), -Coordinates of a second point, localized on the x, y and z-axis, corresponding to the end of the cylindrical diffuser(xout, yout, zout).

  photon history, another photon propagation thread is repeated in the same kernel. Every kernel repeats the entire simulation process until the expected number of photons is reached to obtain an accurate result.The output of the simulation is the Nx.Ny.Nz Absorption matrix [photon_weight/bin], which is normalized to obtain the transport matrix [cm -2 ] and multiplied with the power source [W] to finally obtain the fluence rate [W.cm-2 ].

Figure 3 . 6 : 3 .

 363 Figure 3.6: Pipeline of the GPU Monte-Carlo simulation dedicated to interstitial photodynamic therapy (iPDT) for glioblastoma treatment.

  tissue modelling in light dosimetry experimentations. The concentration of Intralipid controls the scattering coefficient because the absorption coefficient is negligible. Conversely, India ink only affects the absorption coefficient, and not the scattering one. Thus, mixing with the right quantity of Intralipid liquid and India ink in water solution, a phantom with controlled optical properties can be elaborated. Such phantoms have already been extensively studied, and calculations of quantities to be mixed relied on the multi-center study ofSpinelli et al. [296], published in 2014. In this paper, the same optical phantom based on Intralipid and India ink was characterized using nine different methods in nine different research laboratories from six countries. They evaluated the absorption and scattering properties at three different near infrared wavelengths(633, 750 and 830nm). By exploiting results of each laboratory, reference values at these three wavelengths for the intrinsic absorption coefficient of India ink and the intrinsic reduced scattering coefficient of Intralipid concentrated at 20% were determined with an uncertainty of respectively approximately 2% or better, depending on the wavelength considered, and 1%,. Using these values, they proposed a mathematical model to obtain masses of India ink and Intralipid: is the mass of India ink [g], m il is the mass of Intralipid [g], m water is the mass of water [g], m tot is the total mass of the phantom [g], Ԑ a,ink is the intrinsic absorption coefficient of India ink [m -1 ], Ԑ s,il is the intrinsic scattering coefficient of Intralipid [m -1 ], µ a and µ′ s are the desired absorption and reduced scattering coefficients of the phantom.

  specific experimental design was set-up to characterize the fluence rate emitted by different cylindrical diffusors commonly used in PDT treatment within the created phantom. a. Calibration factor The only variable that can be measured in light dosimetry experiment is the light power [W]. To convert this power value Pmeasured [W] into fluence rate value ϕmedium [W.cm -2 ],

Figure 3 . 7 :

 37 Figure 3.7: Optical bench for calibration factor in air estimation.

Figure 3 . 8 :

 38 Figure 3.8: Experimental plan for light power measurements in the phantom.

Figure 3 . 9 :

 39 Figure 3.9: Grid of power measurements positions.

2 cm - 1 - 1 - 1 -

 111 Scattering coefficient: µs = 74.074 cm -Scattering reduced coefficient: µ 𝑠 ′ = 74.074 cm -Anisotropy coefficient: g = 0.733 In each simulation, 10 7 photons were launched isotropically in the surrounding medium. The cylindrical diffuser modeled was placed at the center of the 40x40x80 slab (see figure 3.10).

Figure 3 . 10 :

 310 Figure 3.10: Monte-Carlo simulations of a 40 mm cylindrical diffuser in a homogeneous medium: a) 10 4 photons and b) 10 8 photons were launched during each simulation.

Figure 3 . 11 :

 311 Figure 3.11: Diagram showing the normalized power measured (in nW) according to the incident angle between the probe and the laser beam.

Figure 3 . 12 :

 312 Figure 3.12: Normalized illumination profile of the five cylindrical diffusers described in the table 3.1. The power emitted by the five cylindrical diffusers was measured along each tip.

Figure 3 . 13 :

 313 Figure 3.13: Mean relative deviation according to relative distances to the fiber tip.

  4: K620 GPU features provided by the "Device Query" program sample developed by NVIDIA. a. Optimization of the execution speed Several configurations have been tested to obtain the fastest execution of Monte-Carlo simulation regarding the number of blocks and threads used per block. A homogeneous medium has been implemented with the following optical properties:

Figure 3 . 14 : 1 -

 3141 Figure 3.14: Relative errors of several Monte-Carlo simulations using different total numbers of photons generated (10 5 to 5.10 7 ).

  .............................................................................................................. 140

  may be achievable. This modality lies on the insertion of optical fibers directly into the tumor targeted. Because of strong light absorption in biological tissues, the diffusing sources located at the tip of the optical fibers have to be placed near or inside the tumor under stereotactic conditions to maximize the treated volume. Preliminary results have shown positive outcomes[118]. The treatment strategy (see figure4.1) includes the PDT treatment after its planning with the help of a TPS. Using imagery (MRI and CT-Scan modalities mainly), the GBM is segmented and classified to identify organs at risk and areas to treat. Then, cylindrical diffusers insertion is simulated within the Leksell stereotactic frame of reference. Once the treatment is planned, the procedure starts with the administration of the precursor of the PS (5-ALA) several hours before the illumination in order to have a sufficient concentration of PS (PpIX) to provoke a PDT effect on tumor cells. Optical fibers are then inserted in patient's brain and illuminate surrounding tissues during a given treatment time as planned before the procedure. Once the PDT procedure ended, radiation therapy and chemotherapy are delivered to the patient to decrease recurrence risks (standard of care).

Figure 4 . 1 :

 41 Figure 4.1: Surgical workflow including the 5-ALA iPDT procedure.

- 2 ]

 2 . In this paragraph, a description of the developed software, the validation of the stereotactic registration and dosimetry model implemented is presented [312]. a. iPDT-Plan: a proof of concept This software, named iPDT-Plan, includes stereotactic registration of medical images (digital imaging and communications in medicine, DICOM), manual delineation of the tumor target, light source insertion and dosimetry modeling. The software was implemented with the integrated development environment C++ Builder XE8 (Embarcadero, Austin, TX, USA) and the environment ArtiMED 4 developed at the ONCO-THAI Lab. This environment provides data management, mathematics and image processing tools. Named iPDT-Plan, this first TPS implementation provides an user-interface composed of four main tabs corresponding to the typical planning workflow (see figure 4.2) [102, 103, 308, 312].

Figure 4 . 2 :

 42 Figure 4.2: Workflow of the iPDT-Plan TPS.

  using the two lengths L1,2 and L3,4 as mentioned in figure 4.4. An addition of 40 mm on the zaxis is due to the extra height of the frame mount. Thus, for a given slice n, the Leksell coordinates pLeskell,i of the four points 1, 2, 3 and 4 can be computed as follows:

Figure 4 .

 4 Figure 4.4: a) 3D rendering showing the Leksell frame to the left and right side of the patient's head. b) MRI exams showing the six dots (enclosed in circles) matching the 6 radioopaque rods included inside the Leksell frame.

Figure 4 . 5 :

 45 Figure 4.5: A 2D cross section map (on the left) and 1D line scan plot (on the right) of fluence rate values for a 40 mm length diffuser computed with Monte-Carlo method. At distances below this MFP´, the analytical light propagation models become inaccurate.

  Higher-order Sobol indices quantify the sensitivity of the output Y variance to the interaction of n-input parameters Xi…n. The second-order Sobol index estimates the output Y sensitivity to Xi and Xj inputs parameter (equation 4.13). 𝑆 𝑖𝑗 = V(E[Y|𝑋 𝑖 , 𝑋 𝑗 ]) -𝑉 𝑖 -𝑉 𝑗 V(𝑌) (4.13) The sensitivity package provided by the Comprehensive R Archive Network 5 was implemented in the statistical software R. This library provides different functions for sensitivity analysis computation, including Sobol indices computation by the Monte-Carlo method. The variance Vi can be developed as (equation 4.14):

Figure 4 .

 4 Figure 4.6: a) A phantom model designed to simulate the head of the patient was fixed on a Leksell frame. b) PMMA rod with 13 radio-opaque marbles spaced 10 mm apart.

Figure 4 .

 4 Figure 4.8 represents the mean relative deviation estimated over the six diffuser lengths (i.e.,

Figure 4 .

 4 Figure 4.9 represents the relative deviation where the mean is estimated of all distances d to the source for each diffuser length simulated, summarized in table 4.3. The mean value of relative deviations equals 2.53% (2.09%).

Figure 4 . 9 :

 49 Figure 4.9: Data points are the mean relative deviations of fluence rate values computed by the analytical method, calculated for all six different source lengths over each distance d (0.75 mm to 10 mm).

Figure 4 .

 4 Figure 4.10: First-order Sobol indices (𝑆 µ 𝑎 ) for different distances to the light source. These Sobol indices quantify the contribution to the analytical model variance to the effect of absorption coefficient variation.

Figure 4 .

 4 Figure 4.11: First-order Sobol indices (𝑆 µ 𝑠 ′ ) for different distances from the light source. These Sobol indices quantify the contribution of the analytical model variance to the effect of reduced scattering coefficient variation.

Figure 4 .

 4 Figure 4.12: Second-order Sobol indices (𝑆 µ 𝑎 µ 𝑠 ′ ) for different distances to the light source. These Sobol indices quantify the sensitivity of the analytical model variance to the interaction of reduced scattering and absorption coefficients variations.

Figure 5 . 1 :

 51 Figure 5.1: Surgical workflow including 5-ALA PDT early after FGR achieved to maximize the extent of resection (EOR).

Figure 5 . 2 :

 52 Figure 5.2: Workflow of the intraoperative PDT procedure.

Figure 5 . 3 :

 53 Figure 5.3: (a) Expanded view of the Herloon® device: an inflatable balloon (1) is inserted into the trocar body (2); a sluice gate allows fluid insertion (3) and a flap valve selectively close or open the internal passage of the fiber guide while ensuring its tightness. (b) Assembled view of the first part of the device: the Herloon© balloon commercialized by B-Braun® Aesculap®.

Figure 5 . 4 :

 54 Figure 5.4: (a) Expanded view of the fiber guide: a borosilicate glass tube (2) is inserted into the stainless steel tube (1); a stainless steel luer adaptor (3), a PTFE seal for impermeability (5) and a stainless steel spring tab (5) hold the fiber guide against the trocar. (b) Assembled view of the fiber guide.

Figure 5 . 5 :

 55 Figure 5.5: (a) Complete lighting applicator (part one and two assembled) (b) Perspective view (generated by SolidWorks ® ) of the whole system: the fiber guide and the single-use balloon are inserted into the trocar; the whole device is placed into a surgical cavity.

Figure 5 . 6 :

 56 Figure 5.6: Experimental plan for fluence rate measurements inside calf brains.

Figure 5 . 7 :

 57 Figure 5.7: Locations of the 15 measurement positions around the balloon.

Figure 5 . 8 :

 58 Figure 5.8: PpIX absorption spectrum.

  were tracking by video during the test. Using positions of each point, graphs showing the nominal stress according to the stretch were plotted. Then, from each curve, the Yeoh's model has been used to obtain Young's modulus value. The Yeoh's model is a phenomenological model dedicated to hyper-elastic nonlinear materials with nearly incompressible deformations [346, 347].

Figure 5 . 9 :

 59 Figure 5.9: Mechanical tests: 1) cut out of the tensile specimens, 2) marking of the two points of the tensile specimens 3) mechanical bench used for tensile test: 4) the tensile specimen is fixed between two jaws and 5) stretched until 6) it breaks.

Figure 5 .

 5 Figure 5.11: a) Skull 3D printed b) 3D printed mold filled with gel with a (1:1.10) mixing c) and d) phantom of the brain placed inside its skull 3D printed.

Figure 5 .

 5 Figure 5.12: a) and b) The phantom was placed inside a CT-scan. c) The lighting applicator was inserted into the resection cavity and filled with iodine until the balloon's wall reached the boundaries of the cavity.

Figure 5 . 13 :

 513 Figure 5.13: Light distribution homogeneity assessment.

Figure 5 .

 5 Figure 5.14: (a) Picture of the experimental set-up of temperature measurement: the balloon was filled with 500 mL of diffusing solution and immersed in a hot-water bath adjusted to 37°C. (b) Screenshot of the thermographic camera: the laser emitted a continuous power of 2 W for two hours.

  b. Transfer function A power of 2 W was emitted from the output of the 70 mm cylindrical diffuser. Figure 5.15 shows all fluence rate [mW/cm²] curves measured around the device for the different balloon volumes (40 mL to 500 mL) at different depth from balloon's border [mm].

Figure 5 . 15 :

 515 Figure 5.15: Fluence rate [mW/cm²] curves measured in calf brain tissues at different distances from the balloon's wall [mm] and for different balloon volumes.

Figure 5 . 16 :

 516 Figure 5.16: Fluence [J/cm²] computed at different distances from the balloon's wall [mm] by normalizing fluence curves for each balloon volume at 200J/cm².

Figure 5 . 17 :

 517 Figure 5.17: Transfer function of lighting duration (minutes) necessary to obtain a therapeutic fluence value of 25 J/cm² at a distance of 5 mm from the balloon's wall depending on the volume of diffusing solution injected into the device.

Figure 5 . 19 :

 519 Figure 5.19: Transfer functions computed from ex-vivo measurements and Monte-Carlo model: normalizing a fluence value of 200J/cm² against the balloon's wall.

Figure 5 . 20 :

 520 Figure 5.20: Plot of the ratio of the fluence rate measured in the brain and the irradiance as a function of the depth for different balloon states.

Figure 5 . 21 :

 521 Figure 5.21: Curves obtained from the eleven tensile specimens characterization using extensometer. Black curves represent the data acquisition and color curves are each Yeoh's model fitted to data.

  b. Conformation to the cavity topology Once the phantom placed in the CT-scan, the balloon was filled incrementally with 40, 60, 70, 80 and 90 mL of diluted iodine solution. The volume of 70 mL of iodine solution injected into the device allowed a good fitting of the cavity resection (see figure 5.22).

Figure 5 . 22 :

 522 Figure 5.22: Multiplanar reconstruction of the CT-scan imaging. The balloon was filled with 70 mL of iodine solution.

Figure 5 . 23 :

 523 Figure 5.23: Normalized light output around the balloon's wall according to the diffusing solution injected [mL].

Figure 5 . 24 :

 524 Figure 5.24: Patient test: setting up using the manual user and treatment of 3D printed patient's head.

6 :

 6 Treatment parameters of the four patients enrolled in the INDYGO clinical trial.

Figure 5 . 25 :

 525 Figure 5.25: Operating room during PDT procedure of the third patient treatment.

  ensured strong enthusiasm for developing the Monte-Carlo method for PDT planning [102, 266]. With the assistance of an automatic tissue segmentation method dedicated to GBM [321, 322], the Monte-Carlo method proposed in chapter III can model the fluence rate according to optical coefficients for each voxel of the MRI volume used for planning. Thus, the fluence rate values computed may consider the optical heterogeneity. Although results of the GPU Monte-Carlo validation are in accordance with the AAPM's recommendations, the accuracy of the Monte-Carlo model could be improved by considering the illumination profile of the cylindrical diffusers used in the treatment. Further validation tests with heterogeneous phantoms could enable a deeper validation of the outputs' algorithm.

  device has been characterized using several ex vivo experiments. With these results, a dosimetry model has been generated to easily deduce the treatment duration according indirectly to the size of the cavity surgery. The resulting transfer function displays a linear relation between the quantity of diffusing liquid used to inflate the device in order to fit the cavity's size and the treatment duration. The computed time enables to emit sufficient light to induce a PDT effect at least at 5 mm depth inside surrounding brain tissues of the cavity. This dosimetry model has been validated by Monte-Carlo simulation. The same transfer function has been generated from Monte-Carlo simulation. An error of 9.998% between ex vivo and Monte-Carlo results was obtained.

Figure 1 . 1 : 25 Figure 1 . 2 : 26 Figure 1 . 3 : 28 Figure 1 . 4 : 30 Figure 1 . 5 : 32 Figure 1 . 6 : 34 Figure 1 . 7 : 36 Figure 1 . 8 : 38 Figure 1 . 9 : 39 Figure 1 . 10 : 42 Figure 1 . 11 : 50 Figure 2 . 1 : 54 Figure 2 . 2 : 55 Figure 2 . 3 : 56 Figure 2 . 4 : 57 Figure 2 . 5 : 58 Figure 2 . 6 : 59 Figure 2 . 7 : 61 Figure 2 . 8 : 62 Figure 2 . 9 : 63 Figure 2 . 10 : 64 Figure 2 . 11 : 65 Figure 2 . 12 : 69 Figure 2 . 71 Figure 2 . 14 : 72 Figure 2 . 15 : 73 Figure 3 . 1 : 80 Figure 3 . 2 : 87 Figure 3 . 3 : 93 Figure 3 . 4 : 95 Figure 3 . 5 : 96 Figure 3 . 6 : 99 Figure 3 . 7 :Figure 3 . 8 :Figure 3 . 9 :Figure 3 . 10 :Figure 3 . 11 :Figure 3 . 12 :Figure 3 . 13 :Figure 3 . 14 :Figure 4 . 1 :Figure 4 . 2 :Figure 4 . 3 :Figure 4 . 4 :Figure 4 . 5 :Figure 4 . 6 :Figure 4 . 8 :Figure 4 . 9 :Figure 4 . 10 :Figure 4 . 11 :Figure 4 . 12 :Figure 5 . 1 :Figure 5 . 2 :Figure 5 . 3 :Figure 5 . 4 :Figure 5 . 5 :Figure 5 . 6 :Figure 5 . 7 :Figure 5 . 8 :Figure 5 . 9 :Figure 5 . 10 :Figure 5 .Figure 5 .Figure 5 . 13 :Figure 5 . 15 :Figure 5 . 16 :Figure 5 . 17 :Figure 5 .Figure 5 . 19 :Figure 5 . 20 :Figure 5 . 21 :Figure 5 . 22 :

 112512261328143015321634173618381939110421115021542255235624572558265927612862296321064211652126927121472215733180328733933495359636993738393103113123133144142434445464849410411412515253545556575859510555135155165175519520521522 Figure 1.1: Illustration of the stimulated emission of radiation. .............................................. 25 Figure 1.2: Jablonski diagram of triple state inversion population. ......................................... 26 Figure 1.3: Laser evolution from 80's to current device. ......................................................... 28 Figure 1.4: Jablonski diagram of photophysical processes. ..................................................... 30 Figure 1.5: a) 5-ALA molecular structure, precursor of the photosensitize b) Protoporphyrin IX molecular structure. ............................................................................................................. 32 Figure 1.6: Illustration of the light depth penetration in biological tissue: a) dependency of the wavelength and b) dependency of the size of the light source (taken from [36]). ................... 34 Figure 1.7: Illustration of lighting modalities available for PDT treatments (illustrations taken from [41]). ................................................................................................................................ 36 Figure 1.8: Illustration of cell death mechanism induced by PDT treatment (inspired by [54, 55]). .......................................................................................................................................... 38 Figure 1.9: Gliomas classification according to the world health organization. ...................... 39 Figure 1.10: Decision algorithm of the HGG treatment plan. .................................................. 42 Figure 1.11 : The two main clinical PDT modalities for neurosurgery a) Interstitial PDT: optical fibers are inserted through the skull into the tumor core that contains heterogeneous tissues (necrosis, tumor, and edema). b) Intraoperative PDT: a balloon device is inserted into brain patient to treat borders of the resection cavity. ............................................................... 50 Figure 2.1: Light absorption process by a chromophore of a geometrical cross-section A [m 2 ] and an effective cross-section σa [m 2 ] linked by the absorption efficiency Qa [dimensionless]. .................................................................................................................................................. 54 Figure 2.2: Harmonic oscillator modeling the scattering process. ........................................... 55 Figure 2.3: Relation between wavelength and associated scattering theory. ........................... 56 Figure 2.4: Scattering process by a chromophore. ................................................................... 57 Figure 2.5: Henyey-Greenstein functions for different anisotropic values: plot a) and polar plot b) of the scattering phase function p(θ) [sr -1 ] according to the angle θ (The forward direction along the original photon trajectory is 0° and the backward direction is 180°). ....... 58 Figure 2.6: Mean free path versus reduced mean free path. The term "reduced" enables to extrapolate several progressions of one photon into one mean displacement. ......................... 59 Figure 2.7: Summary of light interactions: reflection, refraction, scattering and absorption in

Figure 5 . 23 :

 523 Figure 5.23: Normalized light output around the balloon's wall according to the diffusing solution injected [mL]. ...........................................................................................................

Figure 5 . 24 :

 524 Figure 5.24: Patient test: setting up using the manual user and treatment of 3D printed patient's head. .........................................................................................................................

Figure 5 . 25 :

 525 Figure 5.25: Operating room during PDT procedure of the third patient treatment. .............

5 .

 5 photons dans la matière, son algorithme est basé sur les lois physiques qui régissent leur propagation. D'un point de vue simplifié, un photon peut être absorbé ou non dans un milieu. S'il ne l'est pas, il est dévié de sa trajectoire initiale. Ces deux notions définissent les deux principales interactions fondamentales de la propagation des photons dans le domaine du visible et sont donc transposées dans l'algorithme de Monte-Carlo. Dans l'ensemble de l'algorithme de Monte-Carlo, des nombres aléatoires sont injectés lors des calculs des interactions. Ces nombres aléatoires permettent de ne pas reproduire le même parcours pour deux photons et ainsi obtenir une simulation pertinente. Chaque photon est défini par plusieurs paramètres tels que sa position, sa direction et son « poids ». Ce poids peut être interprété comme une « barre de vie » où chaque collision diminue sa valeur jusqu'à atteindre la valeur nulle, i.e., l'absorption totale du photon. Comme la plupart des sources lumineuses ont une émission considérée idéalement comme isotropique (c.-à-d., qui émet la lumière de façon homogène, indépendamment de la direction d'émission), la direction initiale est calculée à l'aide d'un nombre aléatoire (équation 3.3). Une fois la direction calculée, une longueur de déplacement est prédite. Si le photon reste dans le même milieu, même après ce déplacement, ce dernier va se déplacer à la position suivante. Dans le cas où le photon change de milieu lors de son déplacement, celui-ci est modifié en fonction des coefficients optiques des deux milieux. Une transmission ou une réflexion du photon peut avoir lieu à l'interface entre les deux milieux. Une fois le photon déplacé, son poids va diminuer en fonction des coefficients optiques du milieu dans lequel il se trouve. Une nouvelle direction va ensuite être calculée. Ces trois parties (« hop-drop-spin » en anglais) vont être répétées jusqu'à ce que le poids du photon soit considéré comme négligeable. Ce n'est qu'alors que le parcours du photon est terminé. Un nouveau photon sera généré jusqu'à atteindre un nombre suffisant de photons simulés pour obtenir une estimation correcte de la quantité de poids déposée dans l'espace. Cette quantité sera ensuite convertie afin d'obtenir l'estimation du débit de fluence (fluence rate en anglais, [W/cm²]) déposé par la source lumineuse.

2 .

 2 d'une boite noire, remplie par le fantôme. Des mesures de puissance autour de chaque fibre ont été réalisées (voir figures 3.8 et 3.9). Ces mêmes fibres ont été modélisées à l'aide du code Monte-Carlo dans un milieu homogène et possédant les mêmes coefficients optiques que ceux du fantôme. L'erreur relative entre les mesures de puissance considérées comme références et les résultats obtenus par simulation Monte-Carlo a été calculée pour chaque mesure. Une erreur relative moyenne de l'ensemble des fibres égale à 14.97 (± 6.16) a été obtenue. Elle se situe dans l'intervalle d'acceptation recommandé par l'AAPM (erreur inférieure à 20%) et permet donc de valider la précision du code Monte-Carlo accéléré par processeur graphique. Concernant les performances du code proposé, une optimisation de l'utilisation des mémoires liées au processeur graphique a été réalisée. Une accélération d'environ 660 a été obtenue en comparaison avec une version similaire exécutée séquentiellement sur processeur (CPU). Cette accélération est cohérente par rapport à celle reportée dans la littérature et permet également des temps de simulation permettant son utilisation en routine clinique. Développement et évaluation d'une preuve de concept d'un système de planification de traitement dédié au traitement photodynamique du glioblastome par voie interstitielle et en conditions stéréotaxiques

3 .

 3 Développement et évaluation d'un nouveau dispositif médical dédié au traitement photodynamique peropératoire du glioblastome La seconde modalité de traitement du GBM, appelée perPDT (intraoperative PDT en anglais), est l'illumination de la cavité opératoire au cours de l'intervention chirurgicale. Cette modalité s'inscrit dans le processus standard de traitement du GBM. Avant l'intervention, le précurseur du PS (l'acide 5-aminolévulinique, 5-ALA) est administré au/à la patient.e dans le cadre de la FGR. En effet, le 5-ALA est actuellement approuvé comme agent de contraste pour la FGR dans le cadre du traitement du GBM chez l'humain en Europe, Asie, Australie et aux Etats-Unis. Une fois le PS (la protoporphyrine IX, PpIX) concentré dans les cellules tumorales (attente de quelques heures), la résection chirurgicale est réalisée, suivie par la FGR. Lorsque la cavité opératoire est illuminée par une lumière bleue (375 -440 nm), le PS émet une fluorescence rougeâtre (635 nm). L'objectif de la FGR est de guider le geste chirurgical afin d'améliorer la qualité de résection. Cependant, bien que macroscopiquement aucun reliquat tumoral ne soit visible en fin de FGR, des cellules tumorales sont toujours

  Germany) sous le nom de Herloon© balloon, est à l'origine un ballon utilisé en chirurgie endoscopique (voir figure 5.3). Il s'agit d'un ballon gonflable et jetable inséré dans un trocart réutilisable. La seconde partie, développée au sein du laboratoire ONCO-THAI, est un guide de fibre lumineux s'insérant dans le trocart et permettant l'insertion et le placement d'une fibre cylindrique diffusante de 70 mm (RD-ML 70, Medlight, Ecublens, Suisse) au centre du ballon (voir figure 5.4). Une fois assemblée, le ballon est gonflé à l'aide d'une solution diffusante (solution intralipidique diluée à 0.1%) permettant une meilleure homogénéisation de la distribution de la lumière autour du ballon (voir figure 5.5). La caractérisation dosimétrique du dispositif a été réalisée sur banc optique afin de simplifier la procédure chirurgicale, notamment l'estimation du temps de traitement. Plus la taille de la cavité est importante, plus le temps d'illumination est important pour déposer la même fluence. Ainsi la durée d'illumination doit être déductible en fonction de la taille de la cavité, et donc de la quantité de solution diffusante injectée dans le ballon. Pour réaliser cette dosimétrie, des mesures expérimentales sur cerveaux de veaux ont été effectuées. Le dispositif était placé dans une boite noire et entouré de cerveaux de veau (voir figure 5.6). Des mesures de puissance ont été réalisées à 15 positions différentes autour du ballon et allant jusqu'à 10 mm de profondeur dans les tissus. Différents volumes de solution diffusante ont été injectés (de 40 mL à 500 mL) permettant de modéliser l'ensemble des volumes de cavité opératoire. Ainsi pour chaque volume, 150 points de mesures ont été effectués. Pour convertir ces mesures de puissance [W] en débit de fluence [w/cm²], le même protocole utilisé pour la validation du modèle Monte-Carlo présenté précédemment a été utilisé. Un facteur de calibration a donc été calculé pour convertir les mesures de puissance dans les cerveaux de veau en débit de fluence. Une fois ces mesures réalisées, un temps de traitement a été calculé afin de déposer une fluence de 200 J/cm 2 correspondant à 25 J/cm 2 pour chaque volume de solution diffusante injecté. Ainsi, une fonction de transfert liant directement la quantité injectée dans le dispositif au temps d'illumination requis pour obtenir un effet PDT jusqu'à 5 mm dans les berges de la cavité a été obtenue. Cette fonction de transfert a également été modélisée via plusieurs modélisations Monte-Carlo. Une erreur moyenne de 9.998% entre les deux fonctions de transfert a été calculée. Plusieurs autres études ont été réalisées pour la caractérisation du dispositif. La conformation du ballon aux bords de la cavité a été modélisée grâce à la création d'un fantôme mécanique. Une séquence IRM peropératoire d'une résection chirurgicale de GBM a été segmentée. Le crâne a été imprimé en polymère AM3300. Le cerveau comportant ainsi une cavité a été formé à partir d'un moule imprimé en 3D. Le gel constituant le cerveau a été étudié par extensométrie afin de reproduire les propriétés mécaniques de tissus cérébraux humains, notamment l'élasticité caractérisée par le module de Young. Une fois créé, le fantôme a été imaginé par scanner avec un dispositif inséré dans la cavité. Avec le volume optimal pour occuper au maximum la cavité opératoire, le dispositif n'est au contact que de 32% de la surface de la cavité. Malgré ce faible taux, l'impact sur le débit de fluence reste modéré.L'impact sur les coefficients optiques de la non-présence du photosensibilisateur (i.e., la PpIX) lors des mesures expérimentales a été évalué. De la même façon, l'homogénéité de la distribution de la lumière autour du ballon, l'étanchéité du dispositif, l'impact thermique du rayonnement laser et la quantité limite de solution diffusante injectable dans le ballon avant déchirure de la paroi ont été étudiés afin de démarrer les procédures pour mettre en place une étude clinique de phase I.Cette étude pilote de faisabilité nommée INDYGO (i.e., intraoperative 5-ALA mediated photodynamic therapy for treating glioblastoma) est actuellement en cours de recrutement. Elle vise à évaluer la faisabilité́ de la réalisation de la PDT peropératoire au cours de la chirurgie d'exérèse du GBM, sans toxicité immédiate inacceptable. Un recrutement de 10 participants est envisagé. Actuellement, cinq patient.e.s ont été traités ; aucun incident technologique ou médical n'a été reporté.ConclusionL'objectif principal de ce travail de thèse a été de fournir les outils technologiques permettant la mise en oeuvre du traitement du GBM par PDT. Afin d'inclure la totalité des patient.es atteint.es du GBM (cas résécables et non-résécables), les deux modalités de traitement (i.e., iPDT et perPDT) ont été étudiées. Bien que nombreuses améliorations des travaux présentés dans ce manuscrit sont d'ores et déjà prévues, ces premières activités de recherche ont permis le transfert des connaissances de la PDT vers la clinique avec les premiers traitements photodynamiques en neurochirurgie en France.

Patients Photosensitizer and posology Follow up status Illumination features Laws

	et al.[91],	5 patients	hematoporphyrin		One GBM and	One frontal diffuser,
	1981		including	derivative, 5 mg/kg	one grade III	400µm, NA: 0.35
			2 GBM	48-72hours prior	died	coupled with argon
				surgery			laser	(630nm),
								300mW during 45
								minutes
	Powers et al.	7 patients	2	mg/kg	of	Recurrence at 2	One	5-25mm
	[115], 1991	including	Photofrin 24 hours	and 27 weeks for	cylindrical diffusor
			2 grade IV	prior to surgery		GBM, and 6 and	for each patient, the
							8 weeks with 2	light dose was about
							stable patients	1000J
	Krishnamurthy	18	2	mg/kg	of	Enhancement of	Interstitial PDT (25
	et al. [116],	gliomas	Photofrin 24 hours	median survival	mm	cylindrical
	2000		including	prior to surgery		of 116.5 days for	diffuser, 6 insertions
			12 GBM				GBM, 493 days	max) with three
							for AA	goups of different
								light dose : 1500J-
								3700J / 3700J-4400J
								/ 4400J-5900J
	Beck et al.	10	20 mg/kg body	Median survival	Insertion of 4-6
	[117,	118],	patients	weight of 5-ALA		was 15 months	cylindrical diffuser
	2007							per patient, Mean
								fluence was 4320-
								11520 J.
	Johansson et al.	5 GBM	20 or			No progression	Insertion of 4-6
	[19], 2013			30 mg/kg body	at 29, 30 and 36	cylindrical diffuser
				weight of 5-ALA		months, death at	per patient
							3 and 9 months
		Table 1.1: Major clinical trials of interstitial PDT treatment of GBM.

I. Light and matter ............................................................................................................. 53 1

  . Optical properties ................................................................................................. 53 2. Radiometric parameters ........................................................................................

61 II. PDT dosimetry: state of the art ..................................................................................... 66 1

  . Predictive dosimetry ............................................................................................. 66 2. Monitoring of the PDT treatment .........................................................................

76 Chapter conclusion ................................................................................................................. 77

  and 𝛻 2 the Laplacien equal to the second derivative in x, y and z (𝛻 2 𝜙(𝑟, 𝑡) =

	𝜕 2 𝜙(𝑟,𝑡) 𝜕𝑥 2 +	𝜕 2 𝜙(𝑟,𝑡) 𝜕𝑦 2 +	𝜕 2 𝜙(𝑟,𝑡) 𝜕𝑧 2

  finite element method represents the most exploited method in PDE solving[153- 159] (see figure2.14). This method discretizes continuous PDE solving into discrete solutions

according to Galerkin's method. A weak formulation of the PDE is solved at each node of the mesh. The Helmholtz equation H (equation 2.23) can be approximated by a linear combination of basis functions 𝜑 multiplied by coefficients 𝛼: 𝐻 ≈ ∑ 𝜑 𝑖 . 𝛼 𝑖 𝑖

(2.29)

  3.16) New direction cosines 𝑢 𝑥,𝑖+1 , 𝑢 𝑦,𝑖+1 , 𝑢 𝑧,𝑖+1 are then computed following calculations:

	𝑢 𝑥,𝑖+1 =	sin 𝜃 . (𝑢 𝑥,𝑖 . 𝑢 𝑧,𝑖 . cos 𝜑 -𝑢 𝑦,𝑖 . sin 𝜑) √1 -𝑢 𝑧,𝑖 2	+ 𝑢 𝑥,𝑖 . 𝑐𝑜𝑠 𝜃
	𝑢 𝑦,𝑖+1 =	sin 𝜃 . (𝑢 𝑦,𝑖 . 𝑢 𝑧,𝑖 . cos 𝜑 + 𝑢 𝑥,𝑖 . sin 𝜑) √1 -𝑢 𝑧,𝑖
		𝑢 𝑥,𝑖+1 = 𝑠𝑖𝑛 𝜃. 𝑐𝑜𝑠 𝜑
		𝑢 𝑦,𝑖+1 = 𝑐𝑜𝑠 𝜃. 𝑠𝑖𝑛 𝜑
	if 𝑢 𝑧,𝑖 ≥ 0	𝑢 𝑧,𝑖+1 = 𝑐𝑜𝑠 𝜃	(3.18)
	else
		𝑢 𝑧,𝑖+1 = -𝑐𝑜𝑠 𝜃
	Using these new direction cosines, a new step size can be generated and photon can be
	considered for a new iteration.

2

+ 𝑢 𝑦,𝑖 . 𝑐𝑜𝑠 𝜃 𝑢 𝑧,𝑖+1 = -sin 𝜃 . cos 𝜑 . √1 -𝑢 𝑧,𝑖 2 + 𝑢 𝑧,𝑖 . cos 𝜃 (3.17)

When 𝑢 𝑥,𝑖 = 0, 𝑢 𝑦,𝑖 = 0, 𝑢 𝑧,𝑖 ± 1, new direction cosines 𝑢 𝑥,𝑖+1 , 𝑢 𝑦,𝑖+1 , 𝑢 𝑧,𝑖+1 are then computed following calculations:

  Flock et al. were the first to introduce the term of hybrid Monte-Carlo method[241].Several previous Monte-Carlo simulations with different optical properties were aggregated in order to extract a function. This function is then injected into a diffusion theory solving to obtain better results.Tinet et al. continued in this idea and proposed a two stages semianalytical Monte Carlo simulation [242, 243]. The first stage was dedicated to the information generator, i.e., to collect total reflectance and transmittance for each scattering event and to reintroduce them into the algorithm to decrease the required number of photons. The second stage used results obtained from the first one to calculate analytically photons propagation. At the same time, Wang et al. proposed a solution closest to the idea of hybridization of the Monte-Carlo method [244, 245]. Monte-Carlo solving was used to simulate the light propagation in areas close to the source and analytical equations were solved at large distance

from the source. Alexandiakis et al. extended this idea to solve the photons migration in the frequency domain [246, 247]. Modulation of phase and amplitude occurred when light was reflected by turbid media and were used to characterize optical properties [248-250]. This discrimination technic was improved in other studies to differentiate high scattering from low scattering regions, isotropic from anisotropic regions or homogeneous from inhomogeneous regions [251-254]. b. CPU parallelism and Internet

  Five different cylindrical diffusers have been used to compare Monte-Carlo to experimental measurements (see table3.1).

	Manufacturer and model	Length of the tip	Diameter	Core	Numerical	Efficiency
	of the fiber	(mm)	(mm)	diameter	aperture	measured
				(µm)	(dimensionless)
	Biolitec, CD 603-20	20	1.65	600	0.37	93%
	Medlight, RD-ML 25	25	1	400	0.37	87%
	Medlight, RD 30	30	0.98	500	0.48	71%
	Medlight, RD-ML 40	40	1	400	0.37	87%
	Medlight, RD-ML 50	50	1	400	0.37	87%

Table 3 .

 3

1: Characteristics of the five cylindrical diffusers used in experiments

These five cylindrical diffusers, from two main optical fibers manufacturers

(Medlight,

Table 3 . 2

 32

	1.055 ± 0.005

: Calibration factor estimation CFmedium for power measurement conversion into fluence rate values.

Table 3 .

 3 3: Mean relative deviation and standard deviation for each cylindrical diffusers.

	13.29 (4.88)

Table 3 .

 3 .4):

	Device 0: "Quadro K620"
	CUDA Driver Version / Runtime Version 7.50 / 4.0
	CUDA Capability Major/Minor version	5.0
	number:
	Total amount of global memory:	2048 MBytes (2147483648 bytes)
	GPU Clock Speed:	1.12 GHz
	Memory Clock rate:	900.00 Mhz
	Total amount of constant memory:	65536 bytes
	Total amount of shared memory per block: 49152 bytes
	Total number of registers available per	65536
	block:
	Maximum sizes of each dimension of a	1024 x 1024 x 64
	block:
	Maximum sizes of each dimension of a	2147483647 x 65535 x 65535
	grid:

  .5:

	32	2026.826 (8.38)	16384
	64	2097.168 (16.25)	32768
	128	2217.3 (17.88)	65536
	256	2345.622 (8.43)	131072
	512	2343.858 (17.29)	262144
	1024	2343.77 (14.80)	524288
	Table 3.5: Mean computing times of one Monte-Carlo simulation (standard deviation),
	according to the number of block used. 512 threads were implemented in each block in all
		simulations.
	Blocks	Mean time (ms)	Threads generated
	1	5364.194 (26.82)	512
	2	2708.338 (16.52)	1024
	4	2715.358 (22.20)	2048
	8	2095.402 (16.51)	4096
	16	2124.522 (23.45)	8192

Table 3 . 6

 36

	Threads per block	Mean time (ms)	Threads generated
	8	6676.506 (35.07)	256
	16	4304.988 (24.61)	512
	32	2646.69 (26.53)	1024
	64	1777.518 (25.38)	2048
	128	1787.228 (46.83)	4096
	256	1761.98 (5.83)	8192
	512	2035.524 (25.94)	16384
	Thus, the fastest and most stable configuration for this GPU configuration is to divide the grid
	into 32 blocks with 256 threads per block which generates 8192 threads for each simulation
	(see table 3.6).
	With this configuration, the time required to perform the same Monte-Carlo simulation
	takes approximately 1341.35s. The resulting acceleration and efficiency, defined previously
	by the Amdhal's law, are approximately 661.7991 and 0.0808 respectively. Using this results,
	the proportion of execution time 𝛼 that benefits from parallelized implementation is
	approximately 0.9986.

: Mean computing times of one Monte-Carlo simulation, and standard deviation associated in parentheses, according to the number of block used. 512 threads were implemented in each block in all simulations.

Table 3 .

 3

7: Relative errors computed for five different total numbers of photons generated during one simulation.

Table 3 .

 3 8: Mean computing times of Monte-Carlo simulation computed by both CUDA and

  𝐼𝑀𝐴.𝑥 -(𝐴 0,0 . 𝑃 𝐿𝑒𝑘𝑠𝑒𝑙𝑙.𝑥 + 𝐴 0,1 . 𝑃 𝐿𝑒𝑘𝑠𝑒𝑙𝑙.𝑦 + 𝐴 0,2 . 𝑃 𝐿𝑒𝑘𝑠𝑒𝑙𝑙.𝑧 ) 𝑡 𝑦 = 𝑃 𝐼𝑀𝐴.𝑦 -(𝐴 1,0 . 𝑃 𝐿𝑒𝑘𝑠𝑒𝑙𝑙.𝑥 + 𝐴 1,1 . 𝑃 𝐿𝑒𝑘𝑠𝑒𝑙𝑙.𝑦 + 𝐴 1,2 . 𝑃 𝐿𝑒𝑘𝑠𝑒𝑙𝑙.𝑧 ) 𝑡 𝑧 = 𝑃 𝐼𝑀𝐴.𝑧 -(𝐴 2,0 . 𝑃 𝐿𝑒𝑘𝑠𝑒𝑙𝑙.𝑥 + 𝐴 2,1 . 𝑃 𝐿𝑒𝑘𝑠𝑒𝑙𝑙.𝑦 + 𝐴 2,2 . 𝑃 𝐿𝑒𝑘𝑠𝑒𝑙𝑙.𝑧 )

	(4.8)

Table 4 .

 4 2: Relative deviation computed at different distances from the source; the mean is computed for all source lengths at a given distance d.

	d (mm)	Mean (SD)	Max deviation
	0.75	0.72% (0.75%)	1.98%
	1	0.03% (0.93%)	2.55%
	2	0.61% (3.48%)	9.36%
	3	1.46% (1.53%)	4.32%
	4	0.49% (2.99%)	7.08%
	5	1.69% (1.96%)	5.51%
	6	0.98% (1.69%)	4.70%
	7	2.89% (3.03%)	7.61%
	8	0.25% (5.67%)	15.38%

for each distance d. Mean relative deviation values are summarized in table 4.2. The mean of all mean relative deviations summarized in table 4.2 equals 1.23% (2.48%).

Figure 4.8: Data points are the mean relative deviations of fluence rate values computed by the analytical model, calculated for each distance d (0.75 mm to 10 mm) over all six different source lengths.

Table 4 .

 4 3: Relative deviation computed for different source lengths; the mean is computed for all distances.

	Source length (mm)	Mean (SD)	Max deviation
	10	4.94% (2.38%)	8.00%
	15	4.70% (4.41%)	15.38%
	20	1.05% (1.23%)	4.24%
	30	1.70% (1.71%)	5.56%
	40	0.91% (1.46%)	5.00%

Table 4 .

 4 5: Optical coefficients properties according to tissue type.

	Tissue	Absorption: µa (cm -1 )	Scattering: µs (cm -1 )	Anisotropy: g	Refractive index: n	Wavelength [nm]	Reference
	GBM cells infiltrated	0.2	160	0.875	1.4	633	[75, 99, 117, 188]
	White matter	1.58 0.8	51 409	0.96 0.84	1.4	633 630	[305] [324]
	Grey matter	2.63 0.2	60.2 90	0.88 0.89	1.4	633 630	[305] [324]
	Necrotic core	0.7	147	0.85	1.4	630	[325]
	Edema	0.2	184	0.875	1.4	633	[326]
	Skull	0.4	350	0.92	1.55	650	[327, 328]
	Standard brain	0.29	127	0.9	1.4	660	[327]

  This bi-component gel (e.g., component A and component B) is initially used to protect and make electrical installation watertight. Adapting the mixing of A and B enables to change the Young's modulus of the gel and thus to obtain a Young's modulus close to human brain.

	Soza et al., 2005 [343]	8.196 -8.863	0.452 -0.461
	Budday et al., 2015 [344]	1.389 -1.895	0.5
	Stewart et al., 2017 [345]	0.17 -16.06 (brain tumors)	0.4 -0.496
	Table 5.1: Young's modulus and Poisson's ratio of human brain.
	However, although those studies show Young's modulus and Poisson's ratio of human
	brain, a complete characterization of mechanical properties involves other features mainly
	because human brain does not belong to elastic domain but to hyper-elastic domain.
	b. Gel phantom
	In this study, the Young's modulus is the only evaluation criterion. A gel phantom was
	produced in order to mimic the mechanical properties of human brain and to enable the
	evaluation of the lighting device. According to the literature, a mean Young's modulus
	(without distinction on white and grey matter) is approximately 4.44 kPa and a mean
	Poisson's ratio is approximately 0.47.
	The selected matter to create such a phantom is Power gel® (Cellpack, Villmergen,
	Switzerland).

  .7).

	Fb (dimensionless)	1.0625 (0.005)
	Fp (dimensionless)	1.03 (0.02)
	Fn (dimensionless)	1.55 (0.05)
	Fi (dimensionless)	1.0237 (0.07)
	CFair (cm -2 )	31721
	CFmedia (cm -2 )	55082.98 (6872.20)
	Table 5.2: Calibration factor estimation for power measurement conversion into fluence rate
	values.

Balloon [mL] Fluence rate at balloon border [mW/cm 2 ] Error [%] Monte-Carlo simulation Ex vivo experiments

	Mean error (%)			9.29%
	50	280.82 (18.47)	305.79 (72.86)	8.89%
	100	192.93 (4.86)	215.94 (27.98)	11.93%
	150	175.45 (21.13)	179.52 (7.02)	2.13%
	200	139.06 (6.54)	115.54 (7.42)	16.91%
	250	118.30 (6.42)	101.36 (14.74)	14.44%
	300	98.23 (20.56)	84.65 (15.15)	13.82
	350	86.06 (14.32)	75.56 (6.76)	12.19
	400	66.28 (14.82)	70.50 (11.66)	6.38
	450	59.53 (22.37)	57.19 (5.33)	3.93
	500	53.10 (11.13)	52.02 (11.63)	2.03

Table 5 .

 5 4: Comparison of fluence rates obtained from Monte-Carlo simulations and ex vivo measurements.

Table 5 .

 5 5: Comparison of fluence rates obtained from Monte-Carlo simulations and ex-vivo measurements.Although Young's modulus does not seem to follow a linear regression, the scale of values matches well with Young's modulus observed in the literature. Because the Young's modulus value for the (1:1.10) mixing (4.92 kPa) is the closest of the one obtained for the brain in the literature (4.44 kPa), the mechanical phantom was chosen to be constituted of a(1:1.10) mixing to best mimic mechanical properties of human brain.

	7.21	5.14	8.31	7.40

  photodynamic therapy.Neurosurgery, 1996. 38(3): p. 552-6; discussion 556-7. 124. Akimoto, J., J. Haraoka, and K. Aizawa, Preliminary clinical report on safety and efficacy of photodynamic therapy using talaporfin sodium for malignant gliomas.

	Photodiagnosis Photodyn Ther, 2012. 9(2): p. 91-9.
	125. Kaye, A.H., G. Morstyn, and D. Brownbill, Adjuvant high-dose photoradiation
	therapy in the treatment of cerebral glioma: a phase 1-2 study. J Neurosurg, 1987.
	67(4): p. 500-5.
	126. Perria, C., et al., Photodynamic therapy of malignant brain tumors: clinical results of,
	difficulties with, questions about, and future prospects for the neurosurgical
	applications. Neurosurgery, 1988. 23(5): p. 557-63.
	127. Kostron, H., A. Obwegeser, and R. Jakober, Photodynamic therapy in neurosurgery: a
	review. J Photochem Photobiol B, 1996. 36(2): p. 157-68.
	128. Rosenthal, M.A., et al., Promising survival in patients with high-grade gliomas
	following therapy with a novel boronated porphyrin. J Clin Neurosci, 2003. 10(4): p.
	425-7.
	133. Muragaki, Y., et al., Phase II clinical study on intraoperative photodynamic therapy
	with talaporfin sodium and semiconductor laser in patients with malignant brain
	tumors. J Neurosurg, 2013. 119(4): p. 845-52.
	134. Henyey, L.G. and J.L. Greenstein, Diffuse radiation in the Galaxy. Astrophysical
	Journal 1941. 93: p. 70-93.
	135. Binzoni, T., et al., The use of the Henyey-Greenstein phase function in Monte Carlo
	simulations in biomedical optics. Phys Med Biol, 2006. 51(17): p. N313-22.

129. Schmidt, M.H., et al., Evaluation of photodynamic therapy near functional brain tissue in patients with recurrent brain tumors. J Neurooncol, 2004. 67(1-2): p. 201-7. 130. Stylli, S.S., et al., Photodynamic therapy of high grade glioma -long term survival. Journal of Clinical Neuroscience, 2005. 12(4): p. 389-98. 131. Muller, P.J. and B.C. Wilson, Photodynamic therapy of brain tumors-A work in progress. Lasers in Surgery and Medicine, 2006. 38(5): p. 384-389. 132. Eljamel, M.S., C. Goodman, and H. Moseley, ALA and Photofrin fluorescence-guided resection and repetitive PDT in glioblastoma multiforme: a single centre Phase III randomised controlled trial. Lasers Med Sci, 2008. 23(4): p. 361-7. 136. Welch, A.J., M.J.C. van Gemert, and W.M. Star, Definitions and Overview of Tissue Optics, in Optical-Thermal Response of Laser-Irradiated Tissue, A.J. Welch and M.J.C.

van Gemert, Editors. 2011, Springer Netherlands: Dordrecht. p. 27-64.

Monte-Carlo simulation executed on a 2.80 GHz Intel Xeon CPU

mapReduce: System and method for efficient large-scale data processing [online], available at: http://www.google.com/patents/US7650331 (accessedSeptember 21, 2017)

"ArtiMED: a medical imaging collaborative platform", ONCO-THAI, http://www.oncothai.fr/index.php?option=com_content&view=article&id=6&Itemid=140&lang=en(July 21, 2017)

Sensitivity: Global Sensitivity Analysis of Model Outputs, CRAN [online] https://CRAN.Rproject.org/package=sensitivity (accessedAugust 28, 2017)

Acknowledgments

I would like to express my sincere gratitude to the members of the jury for giving me the honor of reviewing my work.

converted into 3D Leksell frame coordinates to be used with the stereotactic Leksell System (Elekta, Stockholm, Sweden) or for their programming on a neurosurgical robot. For a given treatment configuration (number, type and location of the cylindrical diffusers, power emitted by light sources), dosimetry modeling and lighting duration estimation are computed to provide a fluence rate map [W/cm²] which is superimposed on the images.

Using the fluence rate matrix [W.cm -2 ], treatment duration is estimated to deposit a fluence value [J.cm -2 ] on a target pre-determined by the user. The larger the target is, the longer the treatment duration to obtain a PDT effect on the treated volume is.

b. Stereotactic frame registration

As previously described, the Leksell stereotactic frame (Elekta®, Stockholm, Sweden) is used to localize point in the 3D Leksell space from imaging (CT or MRI). Translation from imaging space to the Leksell space is achieved through an MR-or CT-indicator box clipped onto the frame during the exams. Radio-opaque rods representing the characteristic "N" of the 1.0 0.3 0.2 0.6 0.7 0.7 0.9 1.4 0.5 0.7 1.1 0.6 0.5 Rod 2

1.8 1.2 1.3 1.3 1.5 1.1 0.6 0.7 0.8 1.2 0.9 1.3 1.1 Rod 3 0.5 0.4 0.8 1.0 0.9 0.7 0.7 0.9 0.5 0.5 0.6 0.8 0.7 Table 4.1: Distances [mm] between Leksell coordinates of the radio-opaque marbles obtained by the TPS and the GammaPlan® (Elekta, Stockholm, Sweden). b. Dosimetry evaluation the absorption coefficient is the sum of the calf brain absorption coefficient µ 𝑎,𝑏𝑟𝑎𝑖𝑛 and the PpIX absorption coefficient µ 𝑎,𝑃𝑝𝐼𝑋 : µ 𝑎,𝑡𝑜𝑡𝑎𝑙 = µ 𝑎,𝑏𝑟𝑎𝑖𝑛 + µ 𝑎,𝑃𝑝𝐼𝑋 (5.3) where µ 𝑎,𝑡𝑜𝑡𝑎𝑙 is the absorption coefficient due to the brain and PS absorption (cm -1 ), µ 𝑎,𝑏𝑟𝑎𝑖𝑛

is the absorption coefficient of the calf brain (cm -1 ), and µ 𝑎,𝑃𝑝𝐼𝑋 is the absorption coefficient of PpIX (cm -1 ).

The PpIX absorption coefficient µa,PpIX can be expressed according to the PpIX concentration CPpIX and its molar extinction coefficient εPpIX at the precise wavelength λ:

µ 𝑎,𝑃𝑝𝐼𝑋,𝜆 = 𝜀 𝑃𝑝𝐼𝑋,𝜆 . 𝐶 𝑃𝑝𝐼𝑋 (5.4) where µ 𝑎,𝑃𝑝𝐼𝑋,𝜆 is the PpIX absorption coefficient at wavelength λ (cm (5.5) where µ 𝑎,𝑃𝑝𝐼𝑋,635 is the PpIX absorption coefficient at a wavelength of 635 nm (cm -1 ), µ 𝑎,𝑃𝑝𝐼𝑋,𝐶𝑅𝐴𝑁,635 is the PpIX absorption coefficient of a sample of unknown concentration at a wavelength of 635 nm (cm -1 ), µ 𝑎,𝑃𝑝𝐼𝑋,𝐶𝑅𝐴𝑁,405 is the PpIX absorption coefficient of a sample of unknown concentration at a wavelength of 405 nm (cm -1 ), 𝜀 𝑃𝑝𝐼𝑋,405 is the molar extinction coefficient of PpIX at a wavelength of 405 nm (L.mol -1 .cm -1 ), and 𝐶 𝑃𝑝𝐼𝑋 is the PpIX concentration (mol.L -1 ).

c. Phantom creation

Once the correct gel mixing ratio was obtained, a mold of a brain including a surgical cavity was created to build-up the final brain phantom.

The shape of the phantom was designed from images of a patient harboring a tumor.

To best mimic the usability of the lighting device within a surgical cavity, MR images acquired intraoperatively on the intraoperative MRI of the neurosurgery department of University Hospital of Lille were used.

Post-processing was achieved on the MITK software (Medical Imaging Interaction From the segmentation results, polygon models have been generated and exported in Polygon File Format in order to be processed in the software SolidWorks® (Dassault and a photodynamic effect is expected to be observed at least within 5 mm of the cavity margin.

Regarding the mechanic assessment of the device, the value of the conformation coefficient could appear as weak: approximately 32% of the cavity surface was in contact with the balloon wall with 70 mL of diffusing liquid injected in the device. The conformation was improved with higher quantity of diffusing liquid injected: a maximum of 90 mL of diffusing liquid was injected during the test. However, this quantity was not relevant as clinical case reproduction because the compression of surrounding tissues was too high and could provoke ischemia.

The space between the balloon's wall and brain tissues, up to 15 mm in the present case, is composed of air, which does not absorb or scatter the light emitted. Thus, the dosimetry model achieved with the assumption of the balloon entirely surrounded by brain tissues still remains acceptable. Nevertheless, this impact should be quantified using several Monte-Carlo modeling of cavities imaged with intraoperative MRI. By obtaining the location and the quantity of diffusing solution injected, a 3D model of the device could be reconstructed inside the cavity. Heterogeneity of the light distribution in the cavity's margins could be modeled. autres. Plusieurs parcours de photons peuvent être calculés en parallèle ; un parcours de photon par unité de calcul. Le temps de calcul est alors drastiquement réduit [27,28].

List of tables

II. Contributions

Afin de couvrir l'ensemble des patient.es (admissibles en chirurgie ou non), les deux