Chafik Arar
email: arar.chafik@gmail.com

Hamoudi Kalla
email: hamoudi.kalla@gmail.com

Salim Kalla
email: salim.kalla@univ-batna.dz

Sonia Bendib
email: bendib.ss@univ-batna.dz

Sabrina

Fault-Tolerant Real-Time Scheduling Algorithm for Energy-Aware Embedded Systems

In this paper, we propose a fault-tolerant scheduling approach that achieves low energy consumption and high reliability efficiency. Our scheduling solution is dedicated to multi-bus heterogeneous architectures, which take as input a given system description and a given fault hypothesis. It is based on active redundancy to mask a fixed number k of failures supported in the system, so that there is no need for detecting and handling such failures. In order to maximize the system's reliability, the replicas of each operation are scheduled on different reliable processors. Our solution can maximize reliability and reduce energy consumption when using active redundancy.

I. INTRODUCTION

Embedded systems invade many sectors of human activity, such as medical applications, transportation, energy production, robotics, and telecommunication. The progresses achieved in electronics and data processing improves the performances of these systems. As a result, the new systems are increasingly small and fast, but also more complex and critical, and thus more sensitive to faults. The presence of some faults in these systems, accidental (design, interaction, . . . ) as well as intentional (human, virus, . . . ), are inevitable. Due to catastrophic consequences (human, ecological, and/or financial disasters) that could involve a fault, these systems must be fault-tolerant [START_REF] Suri | Editorial: Special section on dependable real-time systems[END_REF]. A fault can affect either the hardware or the software of the system. In this paper, we consider only processors faults in distributed architectures with buses. A bus is a multipoint connection characterized by a physical medium that connects all the processors of the architecture. As we are targeting embedded systems with limited resources (for reasons of weight, encumbrance, energy consumption, or price constraints), we investigate only software solutions.

We address hardware fault-tolerant approaches based on scheduling algorithms, to tolerate processors faults in distributed architectures. The approach that we propose is our most recent work for building fault-tolerant distributed embedded real-time systems. Prior results have been published in [START_REF] Assayad | Tradeoff exploration between reliability, power consumption, and execution time for embedded systems -the tsh tricriteria scheduling heuristic[END_REF], [START_REF] Girault | A novel bicriteria scheduling heuristics providing a guaranteed global system failure rate[END_REF]. In this paper, we are interested in approaches based on scheduling algorithms that maximize reliability and reduce energy consumption when using active redundancy to tolerate processors faults.

The paper is organized as follows. Section II describes the system model and states the faults assumptions. Section III presents our approach for providing fault-tolerance. Finally, Section IV concludes the paper.

II. MODELS

The algorithm is modeled as a data-flow graph, called algorithm graph and noted ALG. Each vertex of ALG is an operation and each edge is a data-dependence. A datadependence corresponds to a data transfer between a producer operation and a consumer operation. o1 → o2 means that o1 is a predecessor of o2, and o2 is a successor of o1. The architecture is modeled by a non-directed graph, noted ARC, where each node is a processor, and each edge is a bus. We assume that the architecture is heterogeneous and fully connected. Figure 1(b) is an example of ARC, with four processors P 1 , P 2 , P 3 and P 4 , and two buses Bus 1 and Bus 2 .

We assume only hardware components failures and we assume that the algorithm is correct w.r.t. its specification, i.e., it has been formally validated, for instance with model checking and/or theorem proving tools. We consider only transient processors faults. We assume that at most k processors faults can arise in the system, and that the architecture includes more than k processors.

III. THE PROPOSED APPROACH

The solution that we propose to tolerate processors faults is based on the active redundancy of operations. The advantage of the active redundancy is that the obtained schedule is static; in particular, there is no need for complex on-line re-scheduling of the operations that were executed on a processor when the latter fails; also, it can be proved that the schedule meets a required real-time constraint, both in the absence and in the presence of faults. In many embedded systems, this is mandatory.

To tolerate upto k arbitrary processors faults, each operation o of ALG is actively replicated on k+1 processors of ARC. For example, to tolerate two processors faults, three replicas of each operation of the ALG given in Figure 1(a) are scheduled on different processors (see Figure 3). We assume that all values returned by the k+1 replicas of any operation o of ALG are identical. In order to maximize the reliability R of the schedule, we propose to use the Global System Failure Rate per time unit (GSFR) function that we have proposed in [START_REF] Girault | A novel bicriteria scheduling heuristics providing a guaranteed global system failure rate[END_REF]. The GSFR is the failure rate per time unit of the obtained multiprocessor schedule. In our approach, the replicas of each operation are scheduled on the best processors that minimizes GSFR.

The GSFR of scheduling an operation o i , noted Λ, is computed by the following equation:

Λ(S n ) = -log R(S n )/U (S n ) (1)
where, S n is the static schedule at step n of the algorithm, and U (S n ) is the total utilization of the processors.

A. Voltage, frequency, and energy consumption

The maximum supply voltage is noted V max and the corresponding highest operating frequency is noted f max . For each operation, its execution time Exe assumes that the processor operates at f max and V max . The execution time of the operation placed onto the processor p running at frequency f (taken as a scaling factor) is:

Exe(X, p, f ) = Exe(X, p)/f (2)
Concerning the power consumption, we follow the model of Zhu et al. [START_REF] Zhu | The effects of energy management on reliability in real-time embedded systems[END_REF]. For a single operation placed onto a single processor, the power consumption Pc is:

P c = P s + h(P ind + P d ) (3)
where, P d = Cef V 2 f , Ps is the static power (power to maintain basic circuits and to keep the clock running), h is equal to 1 when the circuit is active and 0 when it is inactive, P ind is the frequency independent active power, P d is the frequency dependent active power, Cef is the switch capacitance, V is the supply voltage, and f is the operating frequency.

For a multiprocessor schedule S, we cannot apply directly equation 3. Instead, we must compute the total energy E(S) consumed by S, and then divide by the schedule length L(S):

P (S) = E(S)/L(S) (4)
In our approach, as the k+1 replicas of each operation are scheduled actively on k+1 distinct processors, the energy consumed by the system is maximal. In order to reduce energy consumption, we propose to execute the k+1 replicas of an operation with different frequencies f. As all the k+1 replicas of each operation may have different end execution time (see Figure 3 for the replicas of A), we choose to align the execution time of all the replicas by changing the frequency f of each replica. For example, to reduce energy consumption of the system of Figure 3, we align the first two replicas A 1 and A 2 with A 3 by changing frequencies as shown in Figure 4. The new frequency of each replica o i of o is depend on the end execution time of the last scheduled replica of o.

IV. CONCLUSION

We have proposed in this paper a solution to tolerate several processors faults in distributed heterogeneous architectures with multiple-bus topology. The proposed solution is based on active redundancy, and on GSFR to maximize reliability and to reduce energy consumption. Currently, we are working to evaluate our approach.

Figure 1 ( 3 Fig. 1 .

 131 Fig. 1. Algorithm and architecture graphs.

Figure 2 1 Fig. 2 .

 212 Figure 2 presents an example of a fault free schedule, where operations are not replicated.

2 Fig. 3 .

 23 Fig. 3. A reliable fault tolerant schedule.

  time

2 Fig. 4 .

 24 Fig. 4. The proposed scheme to reduce energy consumption


Jan Unkelbach

Bjoern H Menze

Ender Konukoglu

Florian Dittmann

Matthieu Le

Nicholas Ayache

Helen Shi

Helen A Shih

Radiotherapy planning for glioblastoma based on a tumor

come   L'archive ouverte pluridisciplinaire

the spatial distribution of tumor cells in normal appearing regions of the brain. The target volume for radiotherapy planning can be defined as an isoline of the simulated tumor cell density.

This paper analyzes the model with respect to implications for target volume definition and identifies its most critical components. A retrospective study involving 10 glioblastoma patients treated at our institution has been performed. To illustrate the main findings of the study, a detailed case study is presented for a glioblastoma located close to the falx. In this situation, the falx represents a boundary for migrating tumor cells, whereas the corpus callosum provides a route for the tumor to spread to the contralateral hemisphere. We further discuss the sensitivity of the model with respect to the input parameters. Correct segmentation of the brain appears to be the most crucial model input. We conclude that the tumor growth model provides a method to account for anisotropic growth patterns of glioma, and may therefore provide a tool to make target delineation more objective and automated.

Introduction

Gliomas are the most common primary brain tumors. For glioblastoma, the most aggressive form, median survival is a little more than one year after diagnosis, despite treatment schedules combining surgery, external beam radiotherapy, and chemotherapy. The clinical experience suggests that improvements in radiotherapy alone will not lead to patient cure. On the other hand, radiotherapy is the single most effective therapy. Hence, the benefit of radiotherapy to prolong survival is undoubted [START_REF] Michael D Walker | An analysis of doseeffect relationship in the radiotherapy of malignant gliomas[END_REF][START_REF] Michael D Walker | Evaluation of BCNU and/or radiotherapy in the treatment of anaplastic gliomas: a cooperative clinical trial[END_REF][START_REF] Kristiansen | Combined modality therapy of operated astrocytomas grade III and IV. Confirmation of the value of postoperative irradiation and lack of potentiation of bleomycin on survival time: a prospective multicenter trial of the Scandinavian Glioblastoma Study Group[END_REF]. Consequently, we hypothesize that better targeting of the radiation can improve the efficacy of radiotherapy and lead to prolonged survival or reduced side effects. Glioblastoma differ from many solid tumors in the sense that they grow infiltratively. Instead of forming a solid tumor mass with a defined boundary, glioblastoma are characterized by a smooth gradient of the tumor cell density. It is well known that tumor cells infiltrate the adjacent brain tissue and can be found several centimeters beyond the enhancing tumor mass that is visible on MRI [START_REF] Kelly | Imaging-based stereotaxic serial biopsies in untreated intracranial glial neoplasms[END_REF][START_REF] Watanabe | Magnetic resonance imaging and histopathology of cerebral gliomas[END_REF][START_REF] Sj Price | Improved delineation of glioma margins and regions of infiltration with the use of diffusion tensor imaging: an image-guided biopsy study[END_REF]. Functional imaging modalities including amino acid Positron Emission Tomography (PET), including FET (Fluoro-Ethyl-Tyrosine) and MET (Methionine) [START_REF] Grosu | L-(methyl-11c) methionine positron emission tomography for target delineation in resected high-grade gliomas before radiotherapy[END_REF][START_REF] Miwa | Discrepancy between lesion distributions on methionine pet and mr images in patients with glioblastoma multiforme: insight from a pet and mr fusion image study[END_REF][START_REF] Niyazi | Fet-pet for malignant glioma treatment planning[END_REF][START_REF] Götz | 18F]FET-PET imaging for treatment and response monitoring of radiation therapy in malignant glioma patients a review[END_REF], have potential to improve the definition of the gross tumor volume. However, these modalities also fail to identify areas of low tumor cell infiltration. In current practice, radiotherapy planning is typically based on the gross tumor volume (GTV) visible on MRI. Many practitioners account for the infiltrative growth by expanding the GTV with a 1-3 centimeter margin to form the clinical target volume (CTV), which is irradiated to a homogeneous dose of 60 Gy. The current treatment planning procedure can potentially be improved by accounting for two growth characteristics of gliomas that are currently not or not consistently incorporated in treatment planning:

• An anisotropic growth of gliomas. MRI data as well as histological analysis after autopsy or resection shows that glioma growth is anisotropic [START_REF] Matsukado | The growth of glioblastoma multiforme (astrocytomas, grades 3 and 4) in neurosurgical practice[END_REF][START_REF] Coons | Anatomy and growth patters of diffuse gliomas[END_REF]. This accounts for the complex shapes of the visible tumor and is due to mainly three aspects:

-Anatomical boundaries: The dura, including its extensions falx cerebri and tentorium cerebelli, represents a boundary for migrating tumor cells. Also, except for rare cases of CSF seeding, gliomas do not infiltrate the ventricles.

-Tumor cells infiltrate gray matter much less than white matter.

-Tumor cells seem to migrate primarily along white matter fiber tracts.

• A spatially varying tumor cell density. Most gliomas lack a defined boundary and tumor cells can be found several centimeters away from the T1 contrast enhancing tumor and outside of the peritumoral edema visible on T2. The tumor is rather characterized by a continuous fall-off of the tumor cell density [START_REF] Kelly | Imaging-based stereotaxic serial biopsies in untreated intracranial glial neoplasms[END_REF][START_REF] Watanabe | Magnetic resonance imaging and histopathology of cerebral gliomas[END_REF].

These macroscopic growth characteristics are partly known from histopathological analysis after autopsy or resection. A comprehensive review can be found in [START_REF] Coons | Anatomy and growth patters of diffuse gliomas[END_REF]. Mathematically, these macroscopic patterns of tumor evolution can be modeled by partial differential equations of reaction-diffusion type. In one of the most popular models, the Fisher-Kolmogorov model, tumor growth is described via two processes: local proliferation of tumor cells and migration into neighboring tissue [START_REF] Harpold | The evolution of mathematical modeling of glioma proliferation and invasion[END_REF]. The anisotropic growth due to anatomical boundaries and reduced gray matter infiltration is considered via a segmentation of the brain into white matter, gray matter and cerebrospinal fluid (CSF). Preferential spread of tumor cells along white matter fibers can be modeled based on a cell diffusion tensor derived from Diffusion Tensor Imaging (DTI) [START_REF] Jbabdi | Simulation of anisotropic growth of low-grade gliomas using diffusion tensor imaging[END_REF][START_REF] Clatz | Realistic simulation of the 3-D growth of brain tumors in MR images coupling diffusion with biomechanical deformation[END_REF]. This leads to a phenomenological tumor growth model that estimates the spatial distribution of tumor cells in the brain tissue in regions that appear normal using current imaging modalities. Incorporating the tumor growth model into radiotherapy planning can be approached in three steps:

• The shape of the target volume is modified as to match an isoline of the simulated tumor cell density while the total volume is kept identical to the conventional target based on an isotropic margin.

• The prescription dose within the target is redistributed based on the varying tumor cell density, in order to deliver less dose to regions of low tumor cell density, and (possibly) boost the regions with high tumor cell density.

• Unavoidable dose outside the target is pushed into regions that are more likely to be infiltrated by tumor cells.

This paper primarily addresses the first aspect, i.e. the use of the tumor growth model for automatic delineation of a radiotherapy target volume, taking into account spatial growth characteristics. An accompanying paper [START_REF] Unkelbach | Radiotherapy planning for glioblastoma based on a tumor growth model: implications for dose redistribution[END_REF] discusses the implications of the model for redistribution of dose, i.e. the second aspect.

Previous work on the Fisher-Kolmogorov glioma model

The Fisher-Kolmogorov model has been found to reproduce a wide range of glioma growth characteristics at a macroscopic scale. For a review regarding the development of the model see for example [START_REF] Harpold | The evolution of mathematical modeling of glioma proliferation and invasion[END_REF]. The mathematical properties of the Fisher-Kolmogorov equations are discussed in detail in [START_REF] Murray | Mathematical biology: I. An introduction[END_REF] and [START_REF] Murray | Mathematical Biology II: Spatial Models and Biomedical Applications[END_REF]. A clinical application involves the personalization of the model to the patient specific anatomy (see e.g. [START_REF] Angelini | Glioma dynamics and computational models: a review of segmentation, registration, and in silico growth algorithms and their clinical applications[END_REF] for a review).

In addition, approaches to estimate of model parameters from imaging data have been proposed [START_REF] Harpold | The evolution of mathematical modeling of glioma proliferation and invasion[END_REF][START_REF] Konukoglu | Image guided personalization of reaction-diffusion type tumor growth models using modified anisotropic eikonal equations[END_REF][START_REF] Bjoern | A generative approach for image-based modeling of tumor growth[END_REF], which are based on estimating the velocity of tumor growth and the steepness of the cell density fall-off. The velocity of tumor growth can be estimated from a sequence of MRI images [START_REF] Konukoglu | Image guided personalization of reaction-diffusion type tumor growth models using modified anisotropic eikonal equations[END_REF]; estimating the steepness of the tumor density fall-off is based on comparing the size of the contrast enhancing tumor core to the size of the peritumoral edema. These observable quantities relate to the generic model parameters proliferation rate and diffusion coefficient. The model has, in particular, been used in life time prediction [START_REF] Swanson | A mathematical modelling tool for predicting survival of individual patients following resection of glioblastoma: a proof of principle[END_REF][START_REF] Murray | Glioblastoma brain tumours: estimating the time from brain tumor initiation and resolution of a patient survival anomaly after similar treatment protocols[END_REF]. In addition, several extensions of the model have been suggested. Those extensions integrate diffusion tensor imaging [START_REF] Jbabdi | Simulation of anisotropic growth of low-grade gliomas using diffusion tensor imaging[END_REF][START_REF] Clatz | Realistic simulation of the 3-D growth of brain tumors in MR images coupling diffusion with biomechanical deformation[END_REF], tumor hypoxia [START_REF] Gu | Applying a patientspecific bio-mathematical model of glioma growth to develop virtual [18f]-fmiso-pet images[END_REF], and the effect of radiotherapy [START_REF] Rockne | Predicting the efficacy of radiotherapy in individual glioblastoma patients in vivo: a mathematical modeling approach[END_REF] and chemotherapy on the evolution of the tumor cell density [START_REF] Tracqui | A mathematical model of glioma growth: the effect of chemotherapy on spatio-temporal growth[END_REF][START_REF] Swanson | Quantifying efficacy of chemotherapy of brain tumors with homogeneous and heterogeneous drug delivery[END_REF]. The idea of utilizing the model for target definition in radiotherapy has been suggested by several authors [START_REF] Cobzas | Tumor invasion margin on a riemannian space of brain fibers[END_REF][START_REF] Konukoglu | Extrapolating glioma invasion margin in brain magnetic resonance images: Suggesting new irradiation margins[END_REF][START_REF] Bondiau | Biocomputing: Numerical simulation of glioblastoma growth and comparison with conventional irradiation margins[END_REF], who propose ways of describing tumor infiltration beyond the visible margins of the tumor. Cobzas [START_REF] Cobzas | Tumor invasion margin on a riemannian space of brain fibers[END_REF] defines a distance metric in the brain based on DTI to describe the anisotropy in tumor growth. Konukoglu [29] and Bondiau [START_REF] Bondiau | Biocomputing: Numerical simulation of glioblastoma growth and comparison with conventional irradiation margins[END_REF] uses a travelling wave approximation of the full reaction diffusion model to estimate the spatial distribution of tumor cells in normal appearing brain tissue.

Contributions and organization of this paper

Even though the Fisher-Kolmogorov model has previously been studied in the mathematical biology and modeling community, an application for target delineation in clinical practice requires further research in multiple directions. Our work uses the previously published Fisher-Kolmogorov model and aims at characterizing the potential of the model for a true clinical application in radiotherapy treatment planning. We aim at developing a semi-automatic contouring tool for gliomas that is practically useful. The contributions of this paper are as follows:

• Systematic review of model assumptions: In section 2 we provide a review of the tumor growth model in order to systematically compile the assumptions made to apply the model to target delineation. This includes 1) highlighting the qualitative features of the model that are relevant for radiotherapy planning, 2) identify the model parameters determining those features, and 3) illustrate how the imaging data is used to personalize the model for an application to the patient at hand.

• Use case characterization: We provide a clear illustration of use cases of the tumor growth model for target delineation. In section 3 we perform a detailed case study for a tumor located in proximity to falx and corpus callosum, a situation where it is difficult to consistently account for anatomical boundaries in a manual delineation. In A, we show results for 4 additional cases with varying tumor location. The presented results illustrate the main findings from a study involving 10 glioblastoma patients previously treated at our institution.

• Sensitivity analysis: We discuss the sensitivity of the simulated tumor cell density with respect to the model inputs (section 3.5). In particular, we highlight the need for accurate brain segmentation into white matter, gray matter and anatomical barriers.

• IMRT planning: We compare radiotherapy plans optimized for intensity-modulated radiotherapy (IMRT) for model based versus manual target delineation (section 4). We illustrate to what degree differences in target delineation are translated into differences in the delivered dose distribution.

Furthermore, for one patient we compare the model predicted tumor cell density to the spatial growth patterns of the recurrent tumor (section 5). Finally, section 6 discusses the scope of the work and the potential of a practical application of the tumor growth model for target delineation.

Tumor growth calculations

The purpose of the tumor growth model amounts to estimating the density of tumor cells in the regions of the brain that appear normal on MRI. It is assumed that the growth of the tumor is described phenomenologically by local proliferation of tumor cells and migration into neighboring tissue.

Patient specific input data

The tumor growth calculations are based on MR images routinely acquired in clinical practice. The images incorporated into the simulation process are T1, T2, T2-FLAIR, and T1 post gadolinium. The T1 post gadolinium image shows the vascularized gross tumor volume, and the T2-FLAIR image shows the surrounding edematous region. For illustrative purpose, these two images are shown in figure 1 for a case discussed in detail in this paper.

Data processing

In the first data processing step, all MR sequences are registered to the image with highest spatial resolution. This is done using rigid registration with 6 degrees of freedom, utilizing the function FLIRT [START_REF] Jenkinson | A global optimization method for robust affine registration of brain images[END_REF] as part of the toolbox FSL [START_REF] Smith | Advances in functional and structural MR image analysis and inplementation as FSL[END_REF][START_REF] Woolrich | Baysian analysis of neuroimaging data in FSL[END_REF]. In the second step, a segmentation of the brain was obtained using the multimodal brain tumor segmentation algorithm published in [START_REF] Menze | A generative model for brain tumor segmentation in multi-modal images[END_REF]. The algorithm is an Expectation-Maximization (EM) based segmentation method, which uses a probabilistic normal tissue atlas as spatial tissue prior. For every voxel, it estimates the posterior probability for three normal tissue classes (white matter, gray matter, and CSF1 ), as well as the lesion outlines on T1 post gadolinium and T2-FLAIR. The result of the EM segmentation is augmented as described in [START_REF] Unkelbach | Glioblastoma growth modeling for radiotherapy target delineation[END_REF] in order to facilitate a reliable identification of falx cerebri and tentorium cerebelli as anatomical barriers. The segmentation result for the patient in figure 1 is shown in figure 2. In the last step, the reference MR image is registered to the radiotherapy planning CT using rigid registration. The transformation matrix is saved and later applied to register the simulated tumor cell density to the planning CT.

Figure 2: Segmentation of the brain into contrast enhancing core (white), peritumoral edema, white matter, gray matter, and CSF (black).

Underlying tumor growth model

where ρ is the proliferation rate which is assumed to be spatially constant, and D(r) is the 3 × 3 diffusion tensor which depends on location r. The first term on the right hand side of equation ( 1) is the diffusion term that models tumor cell migration into neighboring tissue. The second term is a logistic growth term that describes tumor cell proliferation. The tumor cell density c(t, r) takes values between zero and the carrying capacity c max . In this paper, the diffusion tensor is constructed as

D(r) = D w • I r ∈ white matter D g • I r ∈ gray matter ( 2
)
where I is the 3 × 3 identity matrix, and D g and D w are scaling coefficients for gray and white matter, respectively. This construction of the diffusion tensor allows for the modeling of reduced gray matter infiltration via a different diffusion coefficient D g < D w . Anatomical constraints are accounted for by boundary conditions: We assume that tumor cells only spread in brain tissue and do not diffuse into any tissue that is not segmented as gray or white matter. 2

Qualitative features of the solution

The asymptotic solution of the Fisher-Kolmogorov equation is given by a traveling wave solution. At the time of diagnosis, the tumor cell density takes high values close to the carrying capacity c max within the enhancing tumor mass seen on the post-gadolinium T1 image. This is motivated by the observation that resected tissue from this area appears as only tumor tissue in a pathological examination [START_REF] Coons | Anatomy and growth patters of diffuse gliomas[END_REF]. At large distance from the enhancing core, the tumor cell density approaches zero. It can be shown that, within the traveling wave approximation, the following two properties hold:

1. The fall-off of the cell density with distance from the enhancing core is described by a characteristic profile that remains approximately fixed while the tumor grows [START_REF] Murray | Mathematical biology: I. An introduction[END_REF][START_REF] Konukoglu | Extrapolating glioma invasion margin in brain magnetic resonance images: Suggesting new irradiation margins[END_REF]. At some distance, the tumor cell density drops approximately exponentially with distance from the visible tumor mass. For a region comprising white matter, the cell density behaves according to

c(|r|) ∝ exp - |r| λ w (3)
where |r| denotes distance from the visible tumor. The parameter λ w is called the infiltration length and denotes the distance at which the cell density drops by a factor of e. It can be shown that λ w relates to the model parameters via

λ w = D w ρ (4)
Within gray matter the cell density drops with an infiltration length λ g = Dg ρ which is believed to be smaller than the infiltration length in white matter.

2. Over time, the tumor front moves outwards into the brain tissue at a constant velocity. The velocity relates to the model parameters via

v w = 2 D w ρ (5)
In a realistic brain geometry, the propagation of the wave is spatially modulated by anatomical boundaries and by different velocities in gray and white matter.

the water diffusion tensor obtained through DTI imaging. Although this idea has been conceptually introduced [START_REF] Clatz | Realistic simulation of the 3-D growth of brain tumors in MR images coupling diffusion with biomechanical deformation[END_REF][START_REF] Jbabdi | Simulation of anisotropic growth of low-grade gliomas using diffusion tensor imaging[END_REF], it is related to several difficulties that are outside the scope of this paper. Problems related to DTI in this context include: (1) even though histopathology provides anecdotal evidence for preferential spread along white matter fibers, there is no data that quantitatively evaluates this effect;

(2) the DTI signal is corrupted within the edematous region where the increased water content leads to reduced fractional anisotropy (FA) and increased apparent diffusion coefficients (ADC). In this paper, we only consider the effects of anatomical boundaries and reduced infiltration of gray matter, because these effects are widely agreed upon, whereas the use of DTI is more speculative at this stage.

Relating the tumor model to imaging data

In order to apply the model to an individual patient, we have to specify how the tumor cell density resulting from the model relates to the MR images. The contrast enhancing region visible in the T1 post gadolinium MR image is related to a breakdown of the blood-brain-barrier (BBB). Even though this only represents a surrogate for the tumor, it is commonly assumed that contrast enhancing volume corresponds to the vascularized, highly cellular part of the gross tumor. For this paper, we therefore assume that the boundary of the contrast enhancing region corresponds to an isoline of the tumor cell density close to the carrying capacity. It is used as the basis for tumor density calculations is described in subsection 2.6.

In current clinical practice, radiotherapy planning for GBM is mostly based on structural MRI imaging (T1 post gadolinium and T2-FLAIR). However, modern imaging modalities including MR spectroscopy (MRS) and positron emission tomography (PET) have been used for glioblastoma imaging and bear potential to improve the characterization of gross and infiltrative disease. In particular, tracers visualizing amino acid metabolism have been studied [START_REF] Grosu | L-(methyl-11c) methionine positron emission tomography for target delineation in resected high-grade gliomas before radiotherapy[END_REF][START_REF] Miwa | Discrepancy between lesion distributions on methionine pet and mr images in patients with glioblastoma multiforme: insight from a pet and mr fusion image study[END_REF][START_REF] Niyazi | Fet-pet for malignant glioma treatment planning[END_REF][START_REF] Götz | 18F]FET-PET imaging for treatment and response monitoring of radiation therapy in malignant glioma patients a review[END_REF], including FET (Fluoro-Ethyl-Tyrosine) and MET (Methionine). These PET scans are not routinely acquired for glioblastoma patients at our institution. However, if these modalities become established, the tumor growth model could be initialized using a PET based delineation of the tumor instead of using the contrast enhancing region.

Calculation of the tumor cell density

For treatment planning, we want to calculate the tumor cell density in the brain at the time the MR images are taken. Using the Fisher-Kolmogorov equation directly is problematic in that context because the initial condition, which leads to a tumor consistent with the images, is unknown. One approach to circumvent this problem consists of using the traveling wave approximation. It is assumed that, at the time of imaging, the tumor cell density is characterized by its asymptotic solution. In this approach, it is assumed that the tumor cell density at the T1 gadolinium outline corresponds to an isoline of the tumor cell density. For this paper, we assume a value of 70% of the carrying capacity. Starting at the outline, the tumor cell density is extrapolated into the brain tissue, assuming the fall-off profile from the traveling wave approximation, and taking into account the anatomical constraints from the brain tissue segmentation. Mathematically, this is achieved by approximating the partial differential equation in (1) by a static Hamilton-Jacobi equation that is solved using a Fast-Marching method. The details of this approach are described in [START_REF] Konukoglu | Extrapolating glioma invasion margin in brain magnetic resonance images: Suggesting new irradiation margins[END_REF].

Model parameters

The tumor growth model depends on four parameters: The diffusion coefficients in white matter D w and gray matter D g , the proliferation rate ρ, and the carrying capacity c max . We now argue that, for our application in radiotherapy planning, the number of parameters can be reduced to two crucial parameters.

The carrying capacity c max can be considered as a scaling factor for the relative normalized tumor cell density that takes values between zero and one. It is straightforward to see that all results presented in this paper are independent of the value of c max and only depend on the relative tumor cell density. We therefore do not have to consider the value of the carrying capacity c max . We are left with the model parameters D w , D g , and ρ. Within the traveling wave approximation, D w and ρ can be expressed through 1) the velocity of the tumor front in white matter obtained from the product of D w and ρ according to equation ( 5), and 2) the infiltration length in white matter obtained from the ratio of D w and ρ according to equation [START_REF] Kelly | Imaging-based stereotaxic serial biopsies in untreated intracranial glial neoplasms[END_REF]. In this paper, we are only interested in obtaining the distribution of tumor cells at the time when treatment starts. Thus, we do not require the velocity of the tumor front but only the infiltration length λ w = D w /ρ. Since the tumor primarily grows in white matter, we focus on the infiltration length in white matter instead of gray matter. If we express the diffusion coefficient in gray matter via the ratio of white and gray matter coefficient, we are left with two parameters of the tumor growth model:

λ w = Dw ρ :
Infiltration length, describing how fast the tumor cell density drops with distance from the visible tumor volume.

D w /D g : Ratio of diffusion coefficients in white matter and gray matter, parameterizing in parts the anisotropy of tumor growth. The ratio of the infiltration lengths is consequently given by λ w /λ g = D w /D g .

In this paper, we are concerned with the spatial definition of the target volume as an isoline of the tumor cell density. In that context, the exact value of the infiltration length λ w is irrelevant because it does not affect the shape of the isolines. Only the abolute values of the tumor cell density associated with the isolines are determined by λ w3 . Thus, the only relevant model parameter for this paper is the ratio D w /D g , which determines the shape of the isolines of the cell density together with the brain segmentation. The literature consistently suggests that tumor cells infiltrate gray matter much less than white matter 4 This suggests a large value for D w /D g ≫ 1. Most illustrations in this paper were obtained for D w /D g = 100. The most appropriate value is however uncertain and we discuss the impact of uncertainties in D w /D g in section 3.5.

Model based target delineation

In this section we demonstrate the use of the tumor growth model for target definition. In sections 3.1-3.5, we provide a detailed discussion of the patient shown in figure 1 in order to illustrate the main findings of the study. We first compare the model-derived target volumes to the manual delineation used in the clinical treatment plan. We then discuss the impact of anatomical boundaries and reduced gray matter infiltration in more detail, and finally discuss sensitivity to model inputs. We analyzed 10 GBM patients with varying tumor locations previously treated at our institution. In section 3.6, we summarize the main findings for these patients; further details are provided in A.

Target volume definition

Figure 4 shows the clinical target volume (dark green contour) and the boost volume (green contour) drawn manually by the physician, i.e. these contours were the basis for the treatment plan that was actually delivered to the patient. The prescribed dose to the boost volume was 60 Gy, the prescribed dose to the CTV was 46 Gy. The boost volume is defined based on a 2 cm isotropic extension of the contrast enhancing lesion. The CTV is defined via a 1.5 cm expansion of the T2-FLAIR abnormality. Both volumes were subsequently trimmed manually to account for anatomical boundaries (dura, ventricles, falx, and tentorium cerebelli). The dose distribution of the 3D conformal treatment plan delivered to the patient is shown in figure 6d.

Using the tumor growth model, the target volume can be defined as an isoline of the tumor cell density. This is illustrated here for the tumor cell density shown in figure 3a obtained for the parameter value D w /D g = 100. The red and the orange contours in figure 4 show the CTV and the boost volume derived from the tumor growth model, respectively. In this example, the target defining isolines are chosen such that the total enclosed volume is equal to the manually delineated target. In the following two subsections, we discuss the differences between manual and model-derived target volumes in detail.

Figure 4: Comparison of manually defined targets (light/dark green) and model derived targets (orange/red). In the clinical treatment plan, the light green volume was prescribed to 60 Gy, the dark green volume to 46 Gy. The yellow contour shows the abnormality on T2-FLAIR.

Accounting for anatomical boundaries

In the manual delineation of the CTV used in the clinical plan, it is incorporated that the falx represents an anatomical barrier for the migration of tumor cells. Hence, the isotropic target expansion was trimmed manually. In the tumor growth model, the falx is modeled via a layer of CSF and is automatically accounted for through the assumption that tumor cells only migrate within white and gray matter. However, the corpus callosum connects the two hemispheres of the brain. The tumor growth model describes the migration of tumor cells through the corpus callosum (see figure 3 and5a). As a consequence, the target volume based on the growth model is extended into the contralateral hemisphere. Figure 5a shows the tumor cell density overlaid on the coronal T1 gadolinium image. This illustrates the three-dimensional modeling of tumor spread via the model, including areas superior to the corpus callosum. This is not consistently accounted for in the manual CTV. In the manually drawn target volumes, the target is slightly extended into the contralateral hemisphere on the slices that show the corpus callosum, but not on the slices located superiorly and inferiorly (figure 5b). In the model derived target volumes, the target is extended further into the contralateral hemisphere, and the spread of tumor cells in superior-inferior direction beyond the corpus callosum is modeled.

Modeling reduced infiltration of gray matter

It is typically assumed that tumor cells infiltrate gray matter much less than white matter, which is modeled by D w /D g ≫ 1. In the simulated tumor cell density in figure 3, this leads to a sharp drop of the cell density towards the cortex. In most regions, this effect does not alter the global shape of the target volume substantially because the thickness of the cortical layer is only several millimeters. This can be seen in figure 3 in the posterior portion of the tumor. Here, the dura is essentially the boundary of the target volume. Technically, some gray matter surrounding the fissures is outside of the target defining isoline. However, these regions are small and are not relevant for radiotherapy planning as further discussed in section 4. For this patient, the most pronounced differences between the model derived target and the manually drawn target, which can be attributed to reduced gray matter infiltration, arise in the region of the lateral sulcus (left-anterior part of the tumor). In this region, the large amount of gray matter surounding the lateral sulcus represents a soft barrier for migrating tumor cells. Thus, the tumor growth model can be used to identify regions of functioning brain tissue that can be taken out of the target volume.

Difference in target volumes

In order to quantitatively compare the difference of manual and model-derived target volumes, we calculate the Dice coefficient, which is given by the volume where both structures overlap, divided by the average volume of the two structures. For the boost and CTV volumes in figure 4 the Dice coefficients are 0.78 and 0.77, respectively. Thus, three quarters of the manual target is also contained in the model-derived target, but a substantial portion of approximately one quarter is different.

Sensitivity to input parameters

The model input parameters that determine the spatial distribution of tumor cells are the segmentation of the brain (into white matter, gray matter and CSF) as well as the parameter D w /D g .

Brain segmentation: Anatomical boundaries like the ventricles and the falx are modeled via boundary conditions. As a consequence, correct segmentation of the brain is a crucial input to the model. This can be illustrated for the falx: Correct modeling of the anatomical boundary in our approach relies on identifying a layer of CSF along the falx. If due to limitations in the automatic segmentation, the two hemispheres are connected through white or gray matter bridges (outside of the corpus callosum), the falx is not correctly established as a barrier. In this case, the growth model will incorrectly predict interhemispherical cell migration outside of the corpus callosum. This problem can in particular arise if, due to a tumor with mass effect, the brain tissue is compressed and the falx is pushed towards the contralateral hemisphere. A reliable application of the model for target definition therefore requires reliable segmentation of the ventricles and the extensions of the dura, i.e. falx cerebri and tentorium cerebelli. In addition, correct modeling of contralateral spread requires a reliable segmentation of the corpus callosum as white matter. This can, for example, be problematic in the superior region of the corpus callosum when the slice distance in the MRI scan is large. Approaches to customize brain segmentation methods for this application can be found in [START_REF] Unkelbach | Glioblastoma growth modeling for radiotherapy target delineation[END_REF].

Gray matter infiltration parameter: In addition to the segmentation, the ratio of white and gray matter diffusion coefficient D w /D g influences the shape of the isolines of the tumor cell density. For D g = 0, gray matter represents a hard boundary and tumor cells only spread in white matter. For D w /D g = 1, tumor cells spread equally in white and gray matter and the shape of the target is solely influenced by anatomical constraints. In figure 3 the simulated tumor cell density is compared for D w /D g = 10 and D w /D g = 100. For smaller D w /D g , the cell density is more washed out (figure 3b) compared to a larger D w /D g where the tumor cell density follows more closely the white matter structure (figure 3a). It has been discussed above that the cortical gray matter has a thickness of only a few millimeters. As a consequence, varying D w /D g has little impact on the global shape of the target volume. The most significant changes for this patient are around the lateral sulcus.

Interdependence: There is some interdependence regarding the sensitivity against D w /D g versus segmentation errors. For small values of D w /D g , the result may be sensitive to errors in the segmentation of CSF, whereas errors in the segmentation between white and gray matter become irrelevant. For large values of D w /D g , the segmentation of CSF becomes less crucial because falx and tentorium cerebelli are surrounded by a layer of gray matter which models the boundary. However, errors in the segmentation of white matter versus gray matter become more crucial.

Target volume comparison for different tumor locations

We retrospectively analysed 10 patients with varying tumor location. All patients were treated with IMRT to a dose of 60 Gy in 30 fractions. For the clinically delivered treatment plan, the CTV was defined using isotropic expansions of 1-2 cm to the T2-FLAIR hyperintensity. Manual editing was performed to trim the CTV expansion where there were clear anatomical barriers to spread such as the falx or dura. A PTV expansion of 3 mm was added. For all patients, the tumor cell density is calculated for a parameter value of D w /D g = 100 based on the method described in section 2. The model-based target volume is defined as the isoline that encompasses the same total volume as the manually delineated CTV used in the clinically applied plan. The main findings are summarized as follows:

• The modeling of anatomical barriers leads to differences in the target volume for tumors located close to corpus callosum and falx cerebri. The use of the tumor growth model suggests a further expansion of the target into the contralateral hemisphere.

• The modeling of reduced gray matter infiltration leads to differences in the target volume for tumors located in proximity to major sulci. This effect was most pronounced in the region of the lateral sulci where the use of the tumor growth model suggests regions that can be excluded from the target volume.

• Dice coefficients between manual and model derived target volumes were calculated for all 10 patients (table 1). The overlap between manual and model-based target volumes was 79% on average, ranging from 74% to 84%. This indicates that approximately 20% of the manual target volume is not contained in the modelbased volume (and vice versa).

• For all patients IMRT treatment plans were optimized using a homogeneous 60 Gy prescription dose to the CTV and the planning formulation in section 4. Dice coefficients for the 95% isodose lines (57 Gy) were calculated in order to compare the differences in the high dose region between manual and model derived radiation plans (table 1). The Dice coefficient was 0.84 on average, ranging from 0.80 to 0.87, indicating that a substantial part of the differences in target volume also translates into differences in the high dose region.

In the appendix A, we provide details for four representative cases. These four cases span a range of different tumor locations, i.e. GBMs in the parietal lobe, temporal lobe, frontal lobe, and bilateral corpus callosum.

Location

Dice CTV contour Dice 57 Gy isoline 1 Parietal (Fig. 4) 0.77 0.82 2 Temporal/Frontal (Fig. 8) 0.74 0.81 3 Parietal (Fig. 10) 0.81 0.82 4 Temporal (Fig. 11) 0.84 0.87 5 Corpus Callosum (Fig. 12

Model-based target delineation without a reference volume

In figure 4 we have defined the model derived target contours such that the overall volumes of the targets equal the volumes of the targets in the clinical plan. This has been done for the sake of comparison of the shapes of the target volumes. In a practical scenario, the model is used for target definition without the existence of a manually drawn reference volume. In this case, the choice of the target defining isoline can be performed by either specifying the target expansion in white matter, or by specifying the average expansion:

1. In the first scenario, the treatment planner is provided with the simulated tumor cell density together with the delineation of enhancing core and edema as well as the brain segmentation. The treatment planner can then manually pick a target defining isoline by measuring the distance to the enhancing core along major sheets of white matter.

2. In the second case, a temporary reference volume can be created automatically as an isotropic expansion around the enhancing core, intersected with the brain tissue mask. The model-derived target volume can then be defined as the isoline that encompasses the same total volume as the temporary reference volume. In this case, the treatment planner has to specify the average margin added to the enhancing core.

Treatment plan comparison

In this section we illustrate the impact of the model derived target volumes on the dose distribution. A treatment plan for intensity modulated radiotherapy was optimized, using 9 equally spaced coplanar 6MV photon beams. Dose calculation was performed using the Quadrant Infinite Beam (QIB) algorithm in CERR [START_REF] Deasy | CERR: a computational environment for radiotherapy research[END_REF], and the optimization of IMRT treatment plans was performed using our own implementation of the L-BFGS quasi-newton method. More specifically, we minimize the following piece-wise quadratic objective function

f (d) = η w o η N η i∈Vη d i -d max η 2 + (6)
+ w u T N T i∈T (d pres i -d i ) 2 + + w o T N T i∈T (d i -d max T ) 2 + (7) + w o H N H i∈H (d i -d max i ) 2 + ( 8
)
The first term (6) denotes overdose objectives for the organs at risk (OAR). For the maximum doses d max η and weighting factors w o η , we use the values summarized in table 2. The second term [START_REF] Grosu | L-(methyl-11c) methionine positron emission tomography for target delineation in resected high-grade gliomas before radiotherapy[END_REF] represents under and over dose objectives for all target voxels. The prescribed dose d pres i is set to 60 Gy for all voxels i in the boost volume and 46 Gy in the CTV. The maximum dose d max T in all target voxels is set to 60 Gy. For all unclassified normal tissues surrounding the target volume (including skull, brain tissue, ventricles), we use a conformity objective in which d max i decreases linearly with the distance of voxel i from the target volumes. Here, we define a generalization of the conformity objective to the case of an inhomogeneous prescription dose. More specifically, for every voxel i we define

d max i = max d low , max j d pres j -z ij d grad (9)
where z ij denotes the euclidean distance of the voxel i from another voxel j that has a non-zero prescription dose d pres j

. The parameter d grad is specified in Gy per cm and describes the desired dose falloff in healthy tissue; d low is a lower dose threshold below which dose is not penalized. The optimization parameters used are summarized in table 2. Table 2: Objective function parameters used for IMRT optimization Figure 6 shows the dose distributions obtained for two treatment plans. The plan in figure 6a is based on the manually drawn target volumes, whereas the plan in figure 6b is based on the target volumes derived from the growth model. For the sake of comparison, the plan based on model-derived targets was optimized by imposing an additional constraint on the integral dose delivered to the patient, i.e. we minimize objective function [START_REF] Sj Price | Improved delineation of glioma margins and regions of infiltration with the use of diffusion tensor imaging: an image-guided biopsy study[END_REF][START_REF] Grosu | L-(methyl-11c) methionine positron emission tomography for target delineation in resected high-grade gliomas before radiotherapy[END_REF][START_REF] Miwa | Discrepancy between lesion distributions on methionine pet and mr images in patients with glioblastoma multiforme: insight from a pet and mr fusion image study[END_REF] subject to the constraint

i∈V d i ≤ d int (10)
where V denotes the set of all voxels in the patient, and d int denotes the integral dose obtained for the treatment plan based on the manual target volumes. Figure 6c shows the difference of the plans in figures 6a and 6b. It can be seen that in plan 6b more dose is delivered to the contralateral hemisphere. Under the constraint that integral dose is not increased, this is compensated for by delivering less dose to other regions, e.g. near the lateral sulcus.

To quantify the difference in dose distributions in figure 6, we calculate the Dice coefficients for the volumes enclosed by selective isodose lines. The Dice coefficients for the 57 Gy and the 46 Gy isodose line evaluate to 0.84 and 0.83, respectively. In section 3.4 we calculated the Dice coefficient between manual and model-derived boost and CTV volumes, which are 0.78 and 0.77, respectively. It is apparent that differences in target volumes only partly translate into differences in dose distributions. The rugged shape of the model-derived target volume that arises from reduced infiltration of gray matter is mostly not of relevance for radiotherapy planning. Due to the physical characteristics of photon beams it is not possible to spare thin layers of gray matter surrounding small fissures.

5 Comparison to follow-up imaging

Scope

Ideally, we would like to have a way to validate the predictions of the growth model regarding the tumor cell density. Histopathological analysis of resected brain tissue would represent one way to validate the extent of tumor cell infiltration predicted by the model. However, it is currently not realistic to obtain such data for untreated human patients at high spatial resolution. Tissue samples taken from biopsies or resected tumors may provide information at selected locations. Usually this is, however, limited to the gross tumor volume; taking biopsy samples from healthy appearing brain tissue, which would be relevant for target delineation, is not performed. An alternative approach towards collecting evidence for the plausibility of the growth model consists in evaluating the spatial progression patterns of recurrent tumors in follow-up MRI images. The advantage of this approach is that this data is widely available since high grade glioma patients typically receive follow-up MR imaging every 2-3 months. On the downside, model validation based on follow-up imaging has limitations. This includes the difficulty in distinguishing recurrent tumor from other radiation induced changes (e.g. radiation necrosis), and the large variation in response to radiation between patients. In addition, deformations due to surgery, tumor mass effects, and radiation related brain atrophy may occur. Addressing these issues in their complexity is outside the scope of this paper. Nevertheless, we hypothesize that the analysis of progression patterns of recurrent tumors can provide some degree of evidence for the plausibility of the tumor growth model. Due to the variability in treatment and treatment response, this is not the case for all patients -but for individual cases.

Patient characteristics

The patient discussed in sections 3 and 4 only had a biopsy taken prior to radiotherapy without significant tumor resection. Also, the patient did not receive additional cycles of chemotherapy or antiangiogenic drugs after the initial course of chemoradiation with concurrent temozolomide. Figure 7 shows two follow-up T2-FLAIR images for this patient. The follow-up images were registered to the reference image using rigid registration. The image in figure 7a was taken 3.5 months after the diagnostic image (figure 1) and 2 months after completion of radiotherapy. The image in figure 7b was taken 7.5 months after the diagnostic image, which was the last image taken before the patient expired. Even though there is no histopathological confirmation, the radiological appearance and the patient's death about 9 months post therapy suggest progressive disease rather than other radiation induced changes. Thus, for the remainder of this section, we make the assumption that the follow-up imaging shows the progression of the tumor. Furthermore, the tumor recurs/progresses centrally, which is the case for the vast majority of cases given standard of care treatments. Figure 6d shows the dose distribution of the clinically delivered 3D conformal treatment plan.

Comparison of follow-up images and simulated tumor infiltration

By considering the delivered dose distribution in figure 6d, we see that the progressive tumor is located almost entirely within the high dose region receiving approximately 60 Gy. Under the assumption that a radiation dose of 60 Gy kills a fraction of the tumors cells, the shape of the isolines after therapy is not altered within the high dose region. In this case, the Fisher-Kolmogorov model predicts that the progression of the tumor approximately follows the isolines of the tumor cell density simulated pre therapy. This motivates to compare the simulated cell density isolines directly to the follow-up images without explicitly modeling the effect of chemoradiation on the cell density quantitatively.

In figure 7 we overlaid the tumor cell density simulated on the diagnostic image on the follow-up T2-FLAIR images. It is seen that the progression of the tumor post therapy follows the isolines of the simulated tumor cell density relatively closely. In particular, the follow-up images show the progression of the tumor through the corpus callosum into the contralateral hemisphere as predicted by the model. The latter is supported by the T1 post contrast images 3.5 and 7.5 months post therapy, which show contrast enhancement in the corpus callosum (T1 images not shown).

6 Discussion

The problem addressed in this work

In this work, we aim at improving and automating target volume delineation for glioblastoma. We investigate the use of a phenomenological tumor growth model based on a reaction-diffusion equation for the tumor cell density. The model has been studied in the literature over the past years. Here, we aim at bringing this model closer to an application in radiotherapy treatment planning. The problem addressed here originates from the infiltrative growth characteristics of glioma and the fact that tumor cells infiltrate brain tissue beyond the abnormality visible on current imaging modalities. In clinical practice, an approximately 2 centimeter wide margin is added around the T2-FLAIR abnormality. It is however known that glioma growth is not isotropic. Tumor growth patterns are influenced by anatomical boundaries as well as reduced infiltration of gray matter. Hence, adding an isotropic margin will unnecessarily expose normal tissue or miss areas of concern. Anatomical constraints can partly be accounted for in the manual definition of the CTV, e.g. by excluding ventricles and regions outside of the dura. However, a consistent modeling of more complex 3D anatomical features is difficult to achieve manually.

Main findings

We performed a retrospective study involving 10 GBM patients treated previously at our institution. For these patients, we compare manually drawn target volumes used in the clinically delivered treatment plan to model-derived target volumes defined through an isoline of the tumor cell density. One goal of the study is to characterize the best use cases for model based target delineation, i.e. we identify anatomical situations in which the growth model leads to more consistent target volumes compared to manual delineation. As demonstrated in the paper, this may be the case for tumors located in proximity to falx and corpus callosum. In this situation, the falx represents a boundary for migrating cells, but the corpus callosum provides a route for the tumor to spread to the contralateral hemisphere. The tumor growth model represents a tool to automate and objectify target definition for such cases, which is difficult to achieve manually. In addition, it is demonstrated that modeling of reduced infiltration of gray matter may influence target delineation near major sulci, where large accumulations of gray matter surrounding the sulci represent a soft boundary for migrating tumor cells. In particular, this effect is observed for tumors located in proximity to the lateral sulcus. By optimizing IMRT treatment plans based on model-derived and manually drawn target volumes, it is demonstrated that the differences in target delineation may translate into differences in the delivered dose distribution.

Scope of the work

Technical versus clinical goal: The work presented has both a technical goal and a clinical goal. The technical goal of this project is to develop automated delineation of the CTV while accounting for anisotropic growth patters. Currently, anatomical barriers are to some degree taken into account manually. In complex situations, this is a time consuming process and often leads to inconsistent CTVs. Automated GTV to CTV expansion will greatly accelerate treatment planning and make treatment volumes more objective. In addition to the technical goal, more accurate identification of the highest risk regions of microscopic disease may translate into a clinical advantage. The large isotropic CTV expansions currently used include neighboring low risk normal brain tissue, which may receive unnecessary and potentially harmful irradiation. However, by using smaller CTV margins, clinicians risk missing areas of tumor extension, which may lead to earlier locoregional failure. Due to the heterogeneity among GBM patients, the clinical advantage will be difficult to demonstrate, however, the technical goal of automating CTV delineation is itself worthwhile to pursue.

The role of model validation: Common criticisms to the use of tumor growth models for radiotherapy planning are that models are oversimplified, that a validation of the model based on clinical data is possible only to a limited extent, or that model parameters are unknown or uncertain. For this model, the analysis of growth patterns of recurrent tumors may provide some indication for the plausibility of the model for individual patients. However, further research on model validation based on follow-up imaging, including MRI, MR Spectroscopy, and PET is needed. Nevertheless, we argue that an application of the model may be helpful in radiotherapy target volume definition despite these limitations. The model parameterizes features of glioma growth that are commonly agreed upon. The model can be viewed as a component in an automatic GTV to CTV expansion tool, which defines a distance measure in the brain that accounts for the patient specific anatomy. It allows treatment planners to incorporate their belief on glioma growth characteristics in multiple stages: In the first stage, anatomical boundaries can be accounted for, which is achieved by segmentation of the ventricles, falx cerebri and tentorium cerebelli. This is a promising application because, clearly, the falx represents a barrier for migrating tumor cells while the corpus callosum allows for inter-hemispheric spread of disease. In the second stage, reduced gray matter infiltration can be incorporated by choosing D w /D g > 1. If integrated into contouring software used in clinical practice, this tool can suggest an initial target volume to the physician, which optionally can be modified manually.

Conclusions

The Fisher-Kolmogorov equation represents a model for the spatial distribution of infiltrating tumor cells in normal appearing brain tissue. It can be applied in radiotherapy target delineation by defining target volumes as isolines of the tumor cell density. The model can incorporate spatial growth characteristics of infiltrating gliomas: the effect of anatomical barriers and a reduced tumor infiltration in gray matter. Therefore, the tumor growth model can provide the basis for an automatic contouring tool that provides the option to account for widely agreed growth patterns of glioma. The approach appears particularly useful for tumors located close to the falx and the corpus callosum.

The most crucial input to the model is the segmentation of the brain tissue, in particular the reliable segmentation of the anatomical barriers falx cerebri and tentorium cerebelli.

A Target volume comparison for different tumor locations

In section 3.6, we summarized the results for model based target delineation for 10 GBM patients that were analyzed. In this appendix, we provide details on four representative cases, which span a range of different tumor locations, i.e. GBMs in the parietal lobe, temporal lobe, frontal lobe, and bilateral corpus callosum. The Dice coefficients between the manual and model-derived target volumes are summarized in table 1, together with the corresponding Dice coefficient between the 95% isodose lines of the IMRT plans. A.1 Multifocal temporal/frontal lobe case

Figure 8 shows a multifocal GBM case involving the left temporal lobe as well as the frontal lobe. Figure 8a shows the coronal T1 post contrast image, revealing the contrast enhancing tumor mass in the temporal lobe. The axial T1 post contrast image in figure 8b shows the lesion in the frontal lobe. The simulated tumor cell density in figure 9a illustrates several features discussed in sections 3 and 4: The tumor growth model describes the steep fall-off of the tumor cell density in gray matter, leading to differences in the target volume around the lateral sulcus (figure 8b). In addition, modeling tumor cell infiltration through the corpus callosum leads to differences in the target volume in the contralateral frontal lobe. Figure 9 shows the IMRT plan comparison for a homogeneous 60 Gy prescription to the manual CTV (c) and the model-derived CTV (b). The figure shows that the differences in the target volume partly translate into differences in the dose distribution. In particular, the dose difference plot in figure 9d shows a lower dose in the lateral sulcus region for the model-based plan, and a higher dose in the contralateral hemisphere close to the corpus callosum.

A.2 Parietal lobe case

A.3 Temporal lobe case

Figure 11 shows a GBM located in the left temporal lobe. As for the other patients, the ventricles and the falx represent a boundary for infiltrating tumor cells, whereas the corpus callosum allows for contralateral spread. However, in this case, the tumor is located relatively far from the corpus callosum. Therefore, the model based target contour does not differ significantly from the manually drawn CTV in the region around the midsagittal plane (figure 11b). The largest differences between manual and model based CTV occur in the regions near the lateral sulcus, where the tumor growth model describes the sulcus with its surrounding gray matter as a boundary (figures 11a/b). Note that the left lateral sulcus is shifted anteriorly due to the tumor with mass effect.

A.4 Corpus callosum case

Figure 12 shows a GBM case involving the bilateral corpus callosum. The simulated tumor cell density in 12c models the fall-off of the tumor cell density in the cortical gray matter surrounding the sulci. This is most prominent in the proximity of major sulci including the lateral sulcus and occipital lobe. The effect of reduced gray matter infiltration leads to differences of up to approximately 1 cm between manual and model derived CTV. In regions where the CTV extends all the way to the dura (see e.g. right posterior region in figure 12a), the differences are insignificant because of the limited dose gradients that can be achieved with therapeutic photon beams. However, in several regions, differences in the CTV contours translate into dose differences in an IMRT plan. This is visible in figure 12b, where the red contour is further expanded into the major fiber tracts of the corona radiata, but the ventricles and regions in the occipital lobe are spared.

Figure 1 :

 1 Figure 1: (a) post-gadolinium T1 weighted image of a glioblastoma located in the left parietal lobe. (b) T2-FLAIR image of the same tumor, showing the surrounding peritumoral edema. (Note that the right side of the image corresponds to the left side of the brain to follow brain imaging conventions.)

  It is assumed that tumor growth is described by two processes: local proliferation of tumor cells and migration of cells into neighboring brain tissue. Mathematically, this is formalized via the Fisher-Kolmogorov equation, a partial differential equation of reaction-diffusion type for the tumor cell density c(r, t) as a function of location r and time t: ∂ ∂t c(r, t) = ∇ • (D(r)∇c(r, t)) + ρc(r, t) 1 -c(r, t) c max

Figure 3 :

 3 Figure 3: Simulated tumor cell density based on the segmentation in figure 2 for parameters λ w = 4.2 mm and D w /D g = 100 (a) and D w /D g = 10 (b) . The color scale refers to the log 10 cell density normalized to the carrying capacity.

Figure 5 :

 5 Figure 5: (a) Log tumor cell density overlayed on the coronal T1 image showing spread of the tumor in the contralateral hemisphere (λ w = 4.2 mm and D w /D g = 100). (b) Comparison of manual CTV (green) and model based CTV (red) (same contours as in figure 4.

Figure 6 :

 6 Figure 6: Comparison of dose distributions (in units of Gy): (a) for a plan based on manual target volumes; (b) for model-derived target volumes. The difference of the dose distribution is shown in (c).Figure (d) shows the planned dose distribution for the 3D conformal plan the patient was treated with.

  Figure 6: Comparison of dose distributions (in units of Gy): (a) for a plan based on manual target volumes; (b) for model-derived target volumes. The difference of the dose distribution is shown in (c).Figure (d) shows the planned dose distribution for the 3D conformal plan the patient was treated with.

Figure 7 :

 7 Figure 7: T2-FLAIR images acquired after therapy, overlaid with isolines of the tumor cell density simulated on the diagnostic images prior to radiotherapy (figure 3a): (a) 3.5 month after diagnosis with isolines for tumor cell densities of 10 -1 , 10 -2 , and 10 -3 , and (b) 7.5 months after diagnosis with isolines for tumor cell densities of 10 -2 , 10 -3 , and 10 -4 .

Figure 8 :

 8 Figure 8: Multifocal GBM located in the left temporal and frontal lobe: (a) coronal T1 post contrast image; (b) axial T1 post contrast image. The manually drawn CTV is shown in green, the model derived CTV in red.

Figure 9 :

 9 Figure 9: (a) Simulated tumor cell density for the patient shown in figure 8; (b-c) Comparison of IMRT dose distributions (in units of Gy) of the plan based on manual target volume (c) and the model-derived target volume (b). The difference of the dose distributions is shown in (d).

Figure

  Figure 10a shows the T2-FLAIR image of a GBM located in the left parietal lobe. Compared to the case discussed in section 3, the tumor is located more posteriorly and closer to the falx. The simulated tumor cell density for D w /D g = 100 is shown in figure 10b. Similar to the case discussed in section 3, the tumor growth model simulates the infiltration of the contralateral side via the corpus callosum. The comparison of model derived target volume (red) to the manually drawn CTV (green) in figure 10a reveals that the tumor growth model suggests further spread into the contralateral hemisphere.

Figure 10 :

 10 Figure 10: GBM located in the left parietal lobe adjacent to the falx: (a) T2-FLAIR image with manual CTV (green) and model derived CTV (red); (b) simulated tumor cell density.

Figure 11 :

 11 Figure 11: GBM located in the left temporal lobe: (a) T1 post contrast image; (b) T2-FLAIR image; (c) simulated tumor cell density. Note that (b) and (c) represent the same slice whereas (a) is located 7.5 mm inferiorly. The manually drawn CTV is shown in green, the model derived CTV in red.

Figure 12 :

 12 Figure 12: GBM located in the corpus callosum: (a) T1 post contrast image; (b) T2-FLAIR image; (c) simulated tumor cell density. Note that (b) and (c) represent the same slice whereas (a) is located 12.5 mm superiorly. The manually drawn CTV is shown in green, the model derived CTV in red.

Table 1 :

 1 Dice coefficients for the CTV contours and the 57 Gy isoline of an IMRT plan for different tumor locations.

	)	0.76	0.82

In this work, we refer to all brain tissue that is neither white matter, gray matter, nor tumor as CSF.

The model can be extended towards anisotropic tumor cell migration within white matter. This is achieved by replacing the identity matrix I in equation (2) by a tensor that is constructed based on

The value of λ w will be crucial for dose prescription considered in the accompanying paper[START_REF] Unkelbach | Radiotherapy planning for glioblastoma based on a tumor growth model: implications for dose redistribution[END_REF].

This holds for the most common case of astrocytomas, not necessarily for oligodentrogliomas (see[START_REF] Coons | Anatomy and growth patters of diffuse gliomas[END_REF] for a review of glioma growth patterns).

Acknowledgment

The project was supported by the Federal Share of program income earned by Massachusetts General Hospital on C06 CA059267, Proton Therapy Research and Treatment Center. Additional support was provided by the ERC Advanced Grant MedYMA.