Ender Konukoglu
email: ender.konukoglu@sophia.inria.fr

Olivier Clatz

Bjoern H Menze

Marc-André Weber

Bram Stieltjes

Emmanuel Mandonnet

Hervé Delingette

Nicholas Ayache

Image Guided Personalization of Reaction-Diffusion Type Tumor Growth Models Using Modified Anisotropic Eikonal Equations

Reaction-diffusion based tumor growth models have been widely used in the literature for modeling the growth of brain gliomas. Lately, recent models have started integrating medical images in their formulation. Including different tissue types, geometry of the brain and the directions of white matter fiber tracts improved the spatial accuracy of reaction-diffusion models. The adaptation of the general model to the specific patient cases on the other hand has not been studied thoroughly yet. In this work we address this adaptation. We propose a parameter estimation method for reaction-diffusion tumor growth models using time series of medical images. This method estimates the patient specific parameters of the model using the images of the patient taken at successive time instances. The proposed method formulates the evolution of the tumor delineation visible in the images based on the reaction-diffusion dynamics therefore it remains consistent with the information available. We perform thorough analysis of the method using synthetic tumors and show important couplings between parameters of the reaction-diffusion model. We show that several parameters can be uniquely identified in the case of fixing one parameter, namely the proliferation rate of tumor cells. Moreover, regardless of the value the proliferation rate is fixed to, the speed of growth of the tumor can be estimated in terms of the model parameters with accuracy. We also show that using the model-based speed we can simulate the evolution of the tumor for the specific patient case. Finally we apply our method to 2 real cases and show promising preliminary results.

I. INTRODUCTION

Brain tumors that start from glial cells, gliomas, form the major class of primary intracranial cancer, [START_REF] Tovi | Mr imaging in cerebral gliomas analysis of tumour tissue components[END_REF], [START_REF] Price | Diffusion tensor imaging of brain tumours at 3t: A potential tool for assessing white matter tract invasion?[END_REF]. These tumors show a high variability in their malignancy. As some of them remain rather benign, i.e. pilocytic astrocytoma, patients suffering from the most malignant forms, glioblastoma multiforme, have an average life expectancy of 1 year [START_REF] Swanson | Virtual brain tumours (gliomas) enhance the reality of medical imaging and highlight inadequacies of current therapy[END_REF]. During the last 20 years there has been vast amount of research on mathematical descriptions of the growth dynamics of gliomas both at microscopic and macroscopic scales. Cellular interactions, effects of intra-cellular dynamics and microscopic invasion have been studied by mathematical models at the microscopic scale [START_REF] Araujo | A history of the study of solid tumour growth: the contribution of mathematical modelling[END_REF]- [START_REF] Kansal | Simulated brain tumor growth dynamics using a three-dimensional cellular automaton[END_REF]. Large scale dynamics such as the average behavior of the tumor, its spatial evolution and its mass effect on the brain have been addressed by macroscopic models [START_REF] Swanson | Virtual brain tumours (gliomas) enhance the reality of medical imaging and highlight inadequacies of current therapy[END_REF], [START_REF] Cristini | Nonlinear simulation of tumor growth[END_REF], [START_REF] Frieboes | Computer simulation of glioma growth and morphology[END_REF], [START_REF] Sanga | Predictive oncology: A review of multidisciplinary, multiscale in silico modeling linking phenotype, morphology and growth[END_REF]- [START_REF] Mohamed | Finite element modeling of brain tumor mass-effect from 3d medical images[END_REF]. These latter models have included in their formulation different anatomical information that are available in medical images. Different tissue types and white matter fiber directions are used as general parameters. However, personalizing these parameters for each patient case and adapting the generic model to specific patient data has not been thoroughly studied yet.

Personalizing the parameters of a tumor growth model for each patient would clinically be important in two aspects: the parameters and the patient-specific model. The parameters of the model could be used to characterize the tumor and help the diagnosis process by providing the speed of growth or quantifying its morphology. The model combined with the personalized parameters, the "patient-specific" model, would give us the opportunity to simulate the evolution of the specific tumor, adapt the therapy to the patient and predict the further evolution.

For the last 5 years specific attention has been given to one class of macroscopic models, the reaction-diffusion models, in the attempt to link tumor growth models to medical images [START_REF] Swanson | Virtual brain tumours (gliomas) enhance the reality of medical imaging and highlight inadequacies of current therapy[END_REF], [START_REF] Clatz | Realistic simulation of the 3d 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], [START_REF] Cruywagen | The modelling of diffusive tumours[END_REF], [START_REF] Swanson | A quantitative model for differential motility of gliomas in grey and white matter[END_REF]- [START_REF] Mandonnet | Computational modeling of the WHO grade II glioma dynamics: principles and applications to management paradigm[END_REF]. These models describe the evolution of the pathology via proliferation of tumor cells and infiltration into the surrounding tissue. Their formulation consists of reaction-diffusion type partial differential equations (PDEs) with the reaction term (first term in Equation 1) representing the proliferation and the diffusion term representing the infiltration [START_REF] Murray | Mathematical Biology[END_REF]. The system

∂u ∂t = ∇ • (D(x)∇u) + ρu(1 -u) (1)
D∇u • -→ n ∂Ω = 0 (2)
is the general building block of such reaction-diffusion models, where u is the tumor cell density, D is a local diffusion tensor (i.e. symmetric positive definite 3x3 matrix), ρ is the proliferation rate, Ω is the brain domain and ∂Ω represents the boundaries of the brain. Equation 1 describes the temporal evolution of the tumor cell density distribution while Equation 2 represents the no-flux boundary conditions.

A. Previous Works on Reaction-Diffusion Type Models

Tracqui et al. proposed one of the first reaction-diffusion models integrating information coming from medical images [START_REF] Tracqui | A mathematical model of glioma growth: the effect of chemotherapy on spatio-temporal growth[END_REF], namely the geometry of the brain and the size of the tumor. In their model, the tumor cells were assumed to diffuse isotropically with a constant rate in the brain and two reaction terms were included to take into account the proliferation and the cell death. They applied their methodology to 2D slices of CT images. A different simulation was run for each slice. In [START_REF] Cruywagen | The modelling of diffusive tumours[END_REF], Cruywagen et al. built on this idea and proposed to use two populations of tumor cells. Swanson et al. in [START_REF] Swanson | A quantitative model for differential motility of gliomas in grey and white matter[END_REF] included the differential motility of tumor cells in reaction-diffusion models based on the observations of Giese et al. in [START_REF] Giese | Migration of human glioma cells on myelin[END_REF] showing that tumor cells migrate faster on myeling sheaths. This new model formulated the higher motility of tumor cells in the white matter than in the grey matter, using the tissue segmentation coming from anatomical images [START_REF] Swanson | A quantitative model for differential motility of gliomas in grey and white matter[END_REF]. They used a spatially varying isotropic diffusion tensor, which took two different values: one in grey and a much higher one in the white matter. Extending this idea of differential motility, Clatz et al. in [START_REF] Clatz | Realistic simulation of the 3d growth of brain tumors in mr images coupling diffusion with biomechanical deformation[END_REF] and later Jbabdi et al. in [START_REF] Jbabdi | Simulation of anisotropic growth of low-grade gliomas using diffusion tensor imaging[END_REF] have included fiber directions (anisotropy) to introduce the directional preference in the diffusion mechanism of tumor cells. In both studies the models rely on anatomical Magnetic Resonance Images (MRIs) and Diffusion Tensor (DT) MRIs to take into account the geometry and the fiber directions. Clatz et al. also have coupled their diffusive model with a linear elastic mechanical model of the brain to be able to describe the mass effect of the tumor. Recently Hogea et al. have extended the anisotropic model and integrated the observation that proliferating tumor cells push each other. This observation is formulated by an extra advection term added to the base equation given in Equation 1 [START_REF] Hogea | Modeling glioma growth and mass effect in 3d mr images of the brain[END_REF]. In order to explain the mass effect of the tumor they also coupled their diffusive model with a mechanical one using nonlinear elasticity. Besides the research on improving the models themselves, several other works showed potential clinical applications of these models by applying them for different purposes. In [START_REF] Cuadra | Atlas-based segmentation of pathological mr brain images using a model of lesion growth[END_REF], [START_REF] Mohamed | Deformable registration of brain tumor images via a statistical model of tumor-induced deformation[END_REF] the authors applied tumor growth models to register anatomical atlases of healthy subjects onto patient images bearing tumors. In some other works [START_REF] Prastawa | Synthetic ground truth for validation of brain tumor mri segmentation[END_REF], [START_REF] Rexilius | A framework for the generation of realistic brain tumor phantoms and applications[END_REF] growth models have been used to create synthetic images bearing tumors which are then used to evaluate segmentation algorithms.

The reaction-diffusion models provide a general framework allowing the integration of information coming from medical images. Once such an integration is achieved the next step is to adapt the model to specific patients data, in other words to personalize the model. This can be done via estimating the parameters of the general model which best simulates the evolution of the tumor observed in the time series of images. The difficulty in this estimation is due to the sparsity of the available information. The reaction-diffusion models describe the temporal evolution of tumor cell density distributions while, in the images we only observe the evolution of the boundaries of the visible part of the tumor. Therefore, the reaction-diffusion models are not directly applicable in this adaptation. The tumor boundary is assumed to correspond to an "iso-density" contour of the tumor cell density distribution [START_REF] Swanson | Virtual brain tumours (gliomas) enhance the reality of medical imaging and highlight inadequacies of current therapy[END_REF], as shown in Figure 1. In this article we use the terms "tumor boundaries", "tumor delineation" and "tumor front" interchangeably to describe the boundary of the visible part of the tumor in the medical images.

B. Previous Works on Parameter Estimation

The task of estimating parameters from time series of images in the context of tumor growth models is a rather unexplored problem. A first attempt was made by Tracqui et al. in [START_REF] Tracqui | A mathematical model of glioma growth: the effect of chemotherapy on spatio-temporal growth[END_REF] where they optimized the parameters of their model by comparing the area of the tumor observed in CT images at different times and the area of the simulated tumor. The drawback of this approach was to use tumor cell densities which requires an initialization of tumor cell density distribution throughout the brain. Since this is not observable in the images assumptions about the tumor cell density distribution has to be made. Moreover, tissue inhomogeneities and observed directional preference of tumor infiltration is not included. In [START_REF] Hogea | An image-driven parameter estimation problem for a reaction-diffusion glioma growth model with mass effects[END_REF], Hogea et al. propose a PDE based method where they invert their model equations and solve a constrained PDE optimization problem to estimate the parameters. Their work takes into account both the growth model for the glioma and its mass effect on the brain tissue. They propose two ways to optimize for the parameters: using tumor cell distributions and using observed deformation based on landmarks. However, there are certain drawbacks of their method. First of all, their first method uses the knowledge of tumor cell density distribution in the brain. As we have mentioned above this information is not available in medical images. The authors mention a probabilistic method to estimate the density distribution while the details of this method may remain open. Moreover, considering that images mask the low density values, estimating the tumor cell distribution is another open problem. Their second method on the other hand links the observed deformation to the tumor growth parameters. This depends strongly on the assumed coupling between mechanical and diffusion models. The uniqueness of the solution for this method is not addressed. Secondly, although they propose the method for 3D they only provide detailed experiments in 1D without including real images, brain geometry, tissue inhomogeneities or other available anatomical information such as fiber directions.

Recently Swanson et al. in [29] proposed a parameter estimation method for the diffusion process in petri-dish. This method is consistent with the observables in the images as it uses the boundaries of the visible tumor rather than tumor cell densities. Authors have derived analytical approximations for the evolution of the tumor delineation for 2D circular growth. Using these solutions they estimated the diffusion coefficient for the petri-dish experiments. The difficulty in applying this method to medical images is that the analytical solutions derived (hence the method itself) assumes radial symmetric growth. This does not have to be the case for brain tumors. The evolution of the tumor is affected by the brain geometry, different tissues and the fiber structures. Besides this, the existence of a reaction term results in a different evolution than pure diffusion. Therefore, this method need to be modified to take into account the reaction term. In another work, Swanson et al.,[START_REF] Swanson | A mathematical modelling tool for predicting survival of individual patients following resection of glioblastoma: a proof of principle[END_REF] address the parameter estimation problem from a different perspective. In this work they use the asymptotic properties of the reaction-diffusion equations and link them to the information in the images for estimating the parameters. The advantage of this approach, like the previously mentioned work, is that it also uses the tumor delineations visible in the images. However, one of the drawbacks of this method is the mapping of the tumor delineations to spheres with the same volume and using these spheres in the computations. Given tissue inhomogeneity, patient specific geometry of both the brain and the tumor and the differential motility of tumor cells, this may be seen as a very strong assumption. As a result, tumors with different parameters are mapped to the same sphere. The second drawback comes from the fact that authors assume tumor cell density values for the extents of the enhancing regions in the MRIs. This is a reasonable assumption unless the values are used in the estimation process. Unfortunately, exact values and their inter-and intra-patient variability are not known. Setting them to arbitrary values introduces a bias on the estimated parameters. Authors also do not provide an analysis of their method. In our preliminary work [START_REF] Konukoglu | Towards an identification of tumor growth parameters from time series of images[END_REF], we proposed a method to estimate the speed of growth of the tumor, also consistent with the observations in the images. Taking into account the brain geometry, tissue inhomogeneity and fiber directions, the estimated speeds in the white and in the gray matter were given in terms of the model parameters. In order to achieve this, starting from the asymptotic properties of reaction-diffusion models we proposed to use a first order anisotropic Eikonal approximation to describe the evolution of the tumor delineation (i.e. the visible tumor front in the images). Using this approximation we formulated the parameter estimation problem. However, the Eikonal approximation proposed in this previous work was first order and did not include higher order effects (i.e. the curvature of the tumor front and the time dependence of the evolution speed) which influence the values of the parameters. Moreover, the formulation of the parameter estimation problem did not utilize all the information available in the images (i.e. the size of the initial tumor). Finally, the method and the estimated parameters for the given results were not analyzed, which is an essential part of the parameter estimation.

In this work, we propose and analyze a parameter estimation method for reaction-diffusion based tumor growth models using time series of medical images. The method is based on the evolution of the tumor delineation rather than tumor cell densities and in this respect it is consistent with the observations in the images. This evolution is formulated using a modified anisotropic Eikonal model which formulates the motion of the tumor delineation taking into account its curved front and the effect of time on its speed. Unlike the previous works, the method presented in this article takes into account tissue inhomogeneities, fiber structures and the real geometry of both the patient's brain and the tumor while keeping consistent with the image information. We also provide extensive analysis of the method and in general the parameter estimation problem for reaction-diffusion models in the context of glioma modeling. Finally, we show preliminary results of the parameter estimation and the "personalization" of growth models on 2 real cases. In Section II, we explain our method, detail the anisotropic Eikonal approximation we use for describing the temporal evolution of the tumor delineation and formulate the parameter estimation problem. In Section III, we evaluate the performance of our method in retrieving the parameters of the reaction-diffusion growth models. We perform thorough analysis of the estimated parameters and the sensitivity to these parameters. In addition to this, we apply our methodology to real images and show preliminary results. Finally, we conclude by discussing these results, the method and the future work in Section IV.

II. METHOD

The parameter estimation methodology and the choice of the estimated parameters depend on the exact formulation of the underlying reaction-diffusion model. In this work we focus on the formulation proposed in [START_REF] Clatz | Realistic simulation of the 3d 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]. However, due to the similarities of reaction-diffusion models the ideas we present here can be easily adjusted to work for other formulations as well. The model for tumor growth proposed in [START_REF] Clatz | Realistic simulation of the 3d growth of brain tumors in mr images coupling diffusion with biomechanical deformation[END_REF] is formulated by the system given in Equations 1 and 2. The diffusion tensor D, in this model, is constructed based on the observations coming from petri-dish experiments and patient images. Giese et al. in [START_REF] Giese | Migration of human glioma cells on myelin[END_REF] showed that glioma cells move faster on the myelin sheaths, the structural element in the white matter. In addition to this, observations from medical images show that tumor cells follow the fiber tracts. On the grey matter on the other hand, tumor cells move slower and observations do not suggest a preferential diffusion direction. Diffusion information coming from the DTI suggests mostly isotropic (close to isotropic) tensors in the grey matter. As a result, in the formulation given in [START_REF] Clatz | Realistic simulation of the 3d growth of brain tumors in mr images coupling diffusion with biomechanical deformation[END_REF] D is constructed as an anisotropic tensor taking into account two different phenomena: differential motility of tumor cells in different tissues and directional preference of tumor cell diffusion in the white matter. The construction of D, which is obtained from the DT-MRI, is as follows:

D(x) = d g I , x ∈ gray matter d w D water , x ∈ white matter (3)
where tumor cells are assumed to diffuse isotropically in the grey matter with a rate d g and diffuse along the fiber tracts in the white matter proportional to the diffusion tensor of the water molecules D water through a coefficient d w . In this construction D water is obtained from DT-MRI and normalized such that the highest diffusion rate in the brain would be 1. We note that in [START_REF] Clatz | Realistic simulation of the 3d growth of brain tumors in mr images coupling diffusion with biomechanical deformation[END_REF] the authors also couple the evolution of the tumor with its mass effect on the brain but in the present study, as a first step, we focus only on the reaction-diffusion part ignoring the mechanical effect. Once the problem for the growing tumor is solved and understood then the parameter estimation can also take into account the mechanical model.

The reaction-diffusion model given by Equations 1, 2 and 3 describes the temporal evolution of local tumor cell densities. As we have noted before, this creates an inconsistency with the observables in the images making the direct application of these models unsuitable for the parameter estimation problem, see Figure 1. In order to solve the parameter estimation problem we need a formulation consistent with the images. The evolution of the tumor delineation should be the phenomenon that is mathematically described instead of the evolution of the tumor cell densities. In section II-A we detail the construction of a formulation, which captures the same dynamics as the reaction-diffusion model but focuses only on the tumor delineation. Once such a formulation is available then one can optimize the parameters using different error measures and optimization schemes. In section II-B we detail our choice for the error measure and the optimization scheme.

A. An Eikonal Approximation for Reaction-Diffusion Models

The asymptotic properties of the reaction-diffusion equations under certain conditions allow us to construct a traveling time formulation for the tumor delineation. In our previous works we have proposed to use such formulations in the context of tumor growth models [START_REF] Konukoglu | Towards an identification of tumor growth parameters from time series of images[END_REF], [START_REF] Konukoglu | A recursive anisotropic fast marching approach to reaction diffusion equation: Application to tumor growth modeling[END_REF]. Here we build on those ideas and improve our formulation.

Reaction-diffusion equations and their asymptotic properties have been well studied in the literature [START_REF] Aronson | Multidimensional nonlinear diffusion arising in population genetics[END_REF], [START_REF] Ebert | Front propagation into unstable states: universal algebraic convergence towards uniformly translating pulled fronts[END_REF]. These properties have been used for different applications [START_REF] Murray | Mathematical Biology[END_REF], [START_REF] Maini | Traveling wave model to interpret a wound-healing cell migration assay for human peritoneal mesothelial cells[END_REF]- [START_REF] Sermesant | An anisotropic multi-front fast marching method for real-time simulation of cardiac electrophysiology[END_REF]. The most important result for our purposes is the existence of a traveling wave solution in the infinite cylinder and in the case of constant coefficients (spatially and temporally constant diffusion tensor D and ρ). Moreover, any initial condition with compact support converges to this solution in time. The traveling wave solution of Equation 1 has the form u(x, t) = u(x -vt) = u(ξ), where ξ is the moving frame and v is the asymptotic speed of this frame, the wavefront. When this solution is plugged into the reaction-diffusion equation we obtain the ordinary differential equation

n ′ Dn d 2 u dξ 2 + v du dξ + ρu(1 -u) = 0, ( 4
)
where n is the direction of motion in the infinite cylinder and the equation describes the traveling wave solution. This is a constant coefficient nonlinear equation and in order to have admissable solutions, the asymptotic speed v should depend on the diffusion tensor D and ρ [START_REF] Murray | Mathematical Biology[END_REF]. This speed is given by the simple relationship: v = 2 √ ρn ′ Dn. Such a property states that all iso-density contours of u at large times under certain conditions will move with a speed of v. Although this information is very useful it is not complete because the convergence of the observed speed to v is slow, in O(1/t). Following the studies of Ebert et al. we can include the effect of this slow convergence and have a time varying estimate of the speed of the moving frame [START_REF] Ebert | Front propagation into unstable states: universal algebraic convergence towards uniformly translating pulled fronts[END_REF]. As the speed of the moving frame converges to v, the profile of u also converges which implies that until convergence different iso-density contours will move at different speeds. The time varying speed v(t) is the estimate for the origin of the moving frame set at u = 0.5 iso-density contour, and different iso-density contours have slightly different convergence properties. However, the effect of the value of the iso-density contour on v(t) is shown to be O(1/t 2 ) and therefore we ignore it [START_REF] Ebert | Front propagation into unstable states: universal algebraic convergence towards uniformly translating pulled fronts[END_REF].

v(t) = √ n ′ Dn(2 √ ρ -3/(2t √ ρ))
The difference between the observed speed of the moving frame, the asymptotic speed and time varying estimate is shown in Figure 2(a). In Figure 2(b) we show the integrals of these speeds starting from the same initial condition to demonstrate the effect of the convergence on the location of the moving frame (which corresponds to the tumor delineation in the context of our work). At this point we can readily formulate a preliminary traveling time formulation for the tumor delineation using The traveling wave has an asymptotic speed shown in dashed curve however, when we observe the speed of the iso-density contour u = 0.5 in time we notice the low rate of convergence to this speed (see the solid curve). An approximation of the speed of the iso-density contour including the convergence effect yields a closer curve to the the observed on (see point-dashed curve). (b) Starting from the same point the integrals of the speed curves, the distances to the initial point, are shown. Notice that we get a much better approximation when we add the convergence effect. All axis are in non-dimensional coordinates.

v(t) as √ ∇T ′ D∇T = 2 √ ρT 4ρT -3 (5)
T (x) = T 0 ∀x ∈ Γ (6)
where T (x) is the function representing the time when the tumor delineation passes through the point x, Equation 6 represents the Dirichlet type boundary condition stating the initial tumor delineation we start evolving from Γ and T 0 is the time elapsed since the tumor has started diffusing until the acquisition of the first image, see Appendix A for details. We observe that if we do not consider the convergence effect and use the asymptotic speed v then the T dependence of the right hand side in the Equation 5 is removed and we could easily replace T 0 by any value we like in Equation 6. The value of T 0 is not available in the images but it can be regarded as another model parameter to be estimated for. The formulation given in Equations 5 and 6 is valid in the infinite cylinder where the evolution is in one direction (in this case the traveling wave is a plane). We can apply this formulation to more general cases (non-planar cases) in 3D by a local linearization assuming that within a voxel the tumor front is planar and the coefficients are constant [START_REF] Konukoglu | A recursive anisotropic fast marching approach to reaction diffusion equation: Application to tumor growth modeling[END_REF]. Then by starting from the initial tumor delineation and sweeping the domain outwards we construct the solution where fronts at each voxel would be patched together linearly. However, such a generalization does not take into account the effect of curvature. In [START_REF] Keener | Mathematical physiology[END_REF] Keener et al. demonstrate a way to take into account the effect of curvature (where they do not take into account the effect of convergence) for slightly curved surfaces in the case of isotropic diffusion. Following the same principles we can derive the general formulation for anisotropic diffusion including the effect of convergence (see Appendix B). This adds a new term in Equation 5 and we obtain

{ 4ρT -3 2 √ ρT -∇ • D∇T √ ∇T ′ D∇T } √ ∇T ′ D∇T = 1 (7)
where the term ∇ • (D∇T / √ ∇T ′ D∇T ) is the effect of the curvature. In the derivation of this term it is assumed that the surface is slightly curved which requires the effect of curvature to be of a lower order than the term 2 √ ρ (see Appendix B). In order to satisfy this we use a saturation function on this effect whose derivative is 1 near 0 and saturates at ±15% of 2 √ ρ. Adding this to our Equation 7 we obtain the final traveling time formulation which describes the evolution of the tumor front based on the reaction-diffusion formalism:

{ 4ρT -3 2 √ ρT -0.3 √ ρ(1 -e -|κ ef f |/(0.3 √ ρ) )} √ ∇T ′ D∇T = 1 (8)
κ ef f = ∇ • D∇T √ ∇T ′ D∇T (9) T (x) = T 0 ∀x ∈ Γ (10)
where we have chosen to use the exponential form for the saturation function but any other choice would work as well. We chose the value of saturation 15% by comparing the analytical solution of the reaction-diffusion equation with the traveling time formulation given in Equation 8for the spherically symmetric growth case [START_REF] Konukoglu | Modeling Glioma Growth and Personalizing Growth Models in Medical Images[END_REF]. This radial analysis showed that saturation values between 5-20% provided the best fit between these two formulations with not much difference in this range. At this point we also notice that the left hand side of Equation 8 can become negative, especially for low values of T . This is due to the fact that the approximations for the time convergence and curvature effects get worse for lower T values [START_REF] Ebert | Front propagation into unstable states: universal algebraic convergence towards uniformly translating pulled fronts[END_REF] and small tumor sizes [START_REF] Keener | Mathematical physiology[END_REF]. In order to overcome this approximation error, in our scheme we do not let the left hand side become less than

{ 4ρT -3 2 √ ρT -0.3 √ ρ(1 -e -|κ ef f |/(0.3 √ ρ) )} ≥ {0.1 √ ρ}, (11)
which serves as the minimum threshold for the speed of the tumor. As a result of this constraint we are sure to have a growing tumor delineation at all times, consistent with the general reaction-diffusion formulation [START_REF] Clatz | Realistic simulation of the 3d 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], [START_REF] Swanson | A quantitative model for differential motility of gliomas in grey and white matter[END_REF]. Equations 8, 9 and 10 combined with the constraint given by Inequality 11, define the formulation describing the evolution of the tumor delineation in 3D, the traveling time formulation. This formulation is based on the hypothesis that the tumor delineation corresponds to an iso-density contour of the tumor cell density u (the value is not specified) whose evolution is defined by the reaction-diffusion model given in Equations 1 and 2.

The traveling time formulation is a second order partial differential equation (a static Hamilton-Jacobi equation) and numerical solvers for such equations have been well studied in the literature [START_REF] Osher | A level set formulation for the solution of the dirichlet problem for hamilton-jacobi equations[END_REF]- [START_REF] Kao | Fast sweeping methods for static hamilton-jacobi equations[END_REF]. In this work we have chosen to adapt the method we proposed in [START_REF] Konukoglu | A recursive anisotropic fast marching approach to reaction diffusion equation: Application to tumor growth modeling[END_REF] due to its fast computation time and its generality on different geometries. This algorithm starts from the initial delineation and sweeps the domain outwards to compute the traveling time values. We provide the further details of this algorithm in Appendix C. We also note that other methods can also be used to solve the traveling time formulation numerically.

The Eikonal model explained above describes the evolution of the tumor delineation visible in the images and captures the same growth dynamics as the reaction-diffusion models. In Figure 3 we show an example evolution simulated using the Eikonal approximation to show that it captures the same growth dynamics as the reaction-diffusion model given in Equation 1. In the figure, for a synthetic tumor we compare the evolution of the iso-density contour u = 0.4 (value consistent with the one proposed in [START_REF] Tracqui | A mathematical model of glioma growth: the effect of chemotherapy on spatio-temporal growth[END_REF] as the imaging threshold) obtained using the reaction-diffusion model (white contours), which uses the tumor cell density throughout the brain, and the evolution obtained using the traveling time formulation (black contours) starting from the innermost white contour (initial location of the u = 0.4 contour). The Euclidean distance between these two evolutions can be given by the distance between their corresponding contours. The average distance between the black and the white contours for the case given in Figure 3 is 0.78 ± 0.69 mm. Considering the usual resolution of such images (1x1x2.6 mm 3 ) we see how similar these two evolutions are. This similarity demonstrates that in the case of medical images where we cannot directly apply the reaction-diffusion models, the traveling time formulation given by Equations 8, 9 and 10 provides us an alternative formulation based on the same dynamics and which can be directly applied to images.

B. The Parameter Estimation Problem

In the reaction-diffusion model given by Equations 1, 2 and 3 we have three different parameters, d w , d g and ρ. In addition to these, in the previous section, by integrating the convergence characteristics of traveling wave solutions into the traveling time formulation we added another parameter T 0 . This gives us 4 parameters to estimate for: (d w ,d g ,ρ,T 0 ). In this work we optimize these parameters such that the evolution we formulate using the traveling time formulation best fits the real evolution observed in the images, which are taken at different times for the same patient.

In order to formulate the parameter estimation problem we need to define an error measure. In a series of N images taken from the same patient at times t 0 , t 1 ..., t N -1 , we have N snapshots of the tumor delineation, one in every image. We do not pose any constraints on the delineations as they could be extracted using any algorithm automatic, semi-automatic or simply by manual delineation. (We only use the end results of this process therefore, we do not go into detail of the segmentation here but refer the interested reader to some recent articles on the subject [START_REF] Corso | Efficient multilevel brain tumor segmentation with integrated bayesian model classification[END_REF].) For a given parameter set, starting from the first time image we can simulate the evolution of the tumor delineation and compare it with the real evolution observed in these N images. We note that the value of t 0 is not known and regarding the time instances we only know the time differences between acquisitions ∆t 0 = 0, ∆t 1 , ..., ∆t N . Combining these we can define where Γ i is the surface enclosing the visible tumor in the image taken at t i and Γ i is the tumor delineation simulated by the traveling time formulation at t i and dist() is the symmetric distance between two surfaces normalized by the surface area. For two given surfaces Γ A and Γ B we define this distance as

C 1 (d w , d g , ρ, T 0 ) = N -1 1 dist(Γ i , Γ i ) (12)
Γ i = {x|T (dw,dg,ρ,T0) (x) = T 0 + ∆t i } (13) with T (x) = T 0 ∀x ∈ Γ 0
̺ (Γ A , Γ B ) = 1 ♯Γ A ∀x∈ΓA ||x -y(x)|| 2 , y(x) = arg y∈ΓB min ||x -y|| 2 (14) dist (Γ A , Γ B ) = 1 2 [̺(Γ A , Γ B ) + ̺(Γ B , Γ A )], (15)
where ||x -y|| is the Euclidean distance between two points in 3D and ♯Γ A represents the number of points in Γ A . In the formulation given in Equation 12we notice that T 0 is the estimate of t 0 . The estimation of t 0 places the time instances, for which we only know the successive differences, on the convergence curve (like the example given in Figure 2). One information we have not used completely in C 1 is the size of the tumor delineation in the first image Γ 0 . In our experiments we observed that in order to correctly map the time instances on the convergence curve (finding the right T 0 based on the convergence curve given in Figure 2) we need to include this size. The inquiry we make is whether it would have been possible to obtain the delineation Γ 0 at the time T 0 using the traveling time formulation if we had started simulating the evolution the moment the tumor started diffusing, namely at T = 0. The assumption we make here is that the tumor started diffusing from a set of isolated small regions. These small regions actually correspond to the avascular masses that start diffusing and speed up with vascularization. In order to include this in our error measure we run the traveling time formulation backwards in time. The simulation starts from Γ 0 and sweeps the domain within the delineation. We do this by solving the Equation 8 within the visible tumor in the first image. This backward evolution provides us a minimum value of T , T min and the corresponding starting point (or a set of points) x min . We notice that if the parameter set d w , d g , ρ, T 0 is consistent with the size of Γ 0 then T min = T 0 . Therefore, the error we need is a function of |T min -T 0 |. In order to have a measure consistent with C 1 we need to convert this time difference into a spatial distance. For this we use the minimum allowable speed value (see Section II-A) v min = 0.1 ρn ′ max D(x min )n max at the point x min , where n max is the principal eigenvector of D(x min ) providing the highest diffusion rate and the factor 0.1 comes from the minimum threshold for the speed of the tumor explained in Section II-A. Using v min we obtain

C 2 (d w , d g , ρ, T 0 ) = (v min (T min )|T min -T 0 |) 2 (16) C = C 1 + C 2 (17)
Combining C 1 and C 2 we obtain the error criteria C we wish to minimize with respect to the model parameters.

The minimization of C is a multidimensional optimization problem and it can be handled using different methods. One important criterion affecting the choice of the minimization algorithm is that explicit derivatives of C with respect to different parameters are not easily available. Another point is that although the parameters have biologically relevant bounds (such as d w , d g , ρ, T 0 > 0) these constraints are not restrictive. Based on these observations we have chosen to use the unconstrained minimization algorithm proposed by Powell in [START_REF] Powell | Uobyqa: unconstrained optimization by quadratic approximation[END_REF]. The attractive feature of this algorithm is that it does not require derivatives of the objective function. Instead, this function's local quadratic approximations are used in the minimization. The algorithm starts by computing several instances of the objective function, constructs the quadratic approximation using these instances and updates the approximation as the minimization proceeds.

The computation time of the proposed parameter estimation method depends on the size of the tumor and more specifically on how much it has grown. The biggest computational load is running the traveling time formulation several times within the optimization algorithm. Depending on the size of the tumor each run can take up to 2 minutes on a 2.4 GHz Intel Pentium machine with 1Gb of RAM. Consequently, the overall optimization takes on the average 40 to 60 minutes which is a short time range, considering the complexity of the problem. This efficiency is obtained due to the advantage of modeling the evolution of the tumor delineation rather than the cell densities.

III. RESULTS

In the evaluation phase of the parameter estimation method, we test the capabilities of the method for retrieving the real parameters of the tumor growth. We first perform tests with synthetic tumors for which the parameters are known and then we apply the method to real cases and show preliminary results. For the tests with synthetic tumors, we construct a dataset of 180 tumors using the reaction-diffusion model composed of 60 different parameter sets at 3 different locations in the brain. The different parameter sets of the model were constructed using different combinations of d w , d g and ρ values given in the table (in the columns to the right of the parameter name) below motivated by the typical values used in the literature [START_REF] Swanson | Virtual brain tumours (gliomas) enhance the reality of medical imaging and highlight inadequacies of current therapy[END_REF], [START_REF] Clatz | Realistic simulation of the 3d 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]. As can be seen from the values for each parameter the final parameter sets cover a large range of growth speed and anisotropy. Each tumor was initialized in a single voxel and grown in the MR image of a healthy subject with a resolution of 1x1x2.6 mm. The diffusion tensor D was constructed using the DT-MRI of the same subject. In order to create the synthetic images of these tumors, we assumed a simple imaging process where a voxel is visualized as tumoral if the number of tumor cells exceeds 40% of the maximum tumor cell capacity the brain parenchyma can handle [START_REF] Tracqui | A mathematical model of glioma growth: the effect of chemotherapy on spatio-temporal growth[END_REF]. For each tumor, the detection and the first image acquisition is made at the moment when the visible tumor size reaches a maximum diameter of 1.5 cm.

d w [mm 2 /

A. Problem of Non-Uniqueness

In the first set of experiments we tried to estimate all the parameters of the reaction-diffusion model (d w , d g , ρ) and the first image acquisition time T 0 (the time elapsed between the emergence of the tumor and its detection) using the traveling time formulation. In these experiments we observed the non-uniqueness of the solution to this problem caused by the coupling between proliferation and diffusion rates and the sparsity of the information contained in the images. The reaction-diffusion model combined with the imaging process can result in very similar evolutions of the tumor delineation with very different parameters. In Figure 4 we show the evolutions of two different tumors (green and red) for which the diffusion and proliferation parameters are given in the accompanying table. The contours with the same color are the delineations of the same tumor in different images taken at successive time instances. The inner contour is the delineation in the first image and the other contours as we go outwards are from the images taken at 200, 300 and 400 days after the time of the first image acquisition respectively. We observe that although the parameters are different the evolutions are almost the same. Quantitatively, the difference between these two evolutions measured by the error criterion C (see Equation 17) is 0.644 mm 2 . On the other hand, the closest tumor delineation evolutions we can get to these ones using the traveling time formulation with the optimum parameters have errors of C = 1.28 mm 2 for the red and C = 1.29 mm 2 for the green tumor. This shows us that with the current resolution of medical images we cannot distinguish between these two parameter sets if we observe either of the evolutions. Therefore, we leave aside the question of estimating separately the diffusion and the proliferation rate.

One observation on the values of the parameters is that between the two cases in Figure 4 the products 2 √ ρd w and 2 ρd g remain almost the same, around 0.1 mm/day and 0.03 mm/day respectively. We have seen in Section II-A that the asymptotic speed of the traveling wave solutions of reaction-diffusion models is given by 2 √ ρn ′ Dn. Therefore, the example shown in the Figure 4 suggests us that the similar volume evolution of these two tumors can be captured and quantified by the asymptotic speed of the model in the white and the grey matter which are both given as functions of the model parameters.

B. Fixing ρ and the 3 Parameter Case

Since estimating all the parameters of the reaction-diffusion equation yielded a non-unique solution (under the given constraints) we turn our attention to the case where we can fix a parameter. The proliferation rate ρ is a microscopic parameter and its coupling with the diffusion rate creates the non-uniqueness of the solution. Here we assume that the value of ρ can be estimated using biopsy results and microscopic analysis. More specifically, we rely on the works showing that the mitotic We observe that although the diffusion and proliferation rates of these tumors are different the evolutions are almost the same. The difference between these evolutions measured using C is 0.644mm 2 which is lower than the minimum error we find by estimating the parameters using the traveling time formulation (C = 1.28 for red and C = 1.29 for green). This shows that we cannot distinguish between these two parameter sets if we observe either of the evolutions. We also observe that the products dwρ and dg ρ are very close for the two tumors. This tells us that although distinguishing between dw, dg and ρ is not obvious estimating the product of these values can be possible.

index (MI) can be computed through the labeling index (LI) and can be linked to the proliferation rate ρ by assuming an average cell cycle duration for the patient [START_REF] Schröder | The relationship between Ki-67 labeling and mitotic index in gliomas and meningiomas: demonstration of the variability of the intermitotic cycle time[END_REF], [START_REF] Johannessen | The clinical value of Ki-67/MIB-1 labeling index in human astrocytomas[END_REF]. We continue our analysis under this assumption. In the first analysis we assume that we know the real value of ρ and fix it in the parameter estimation, once it is fixed the problem becomes uniquely solvable. In this case we are left with three parameters to estimate (d w , d g ) and T 0 .

For each of the synthetic tumors previously described we create a dataset of 3 images, the first image taken at the time of detection and two other images taken at 200 and 400 days after the detection. Using these images and the time difference between acquisitions we estimate the diffusion parameters and T 0 . We show and discuss the obtained estimates based on two different analyses. The first one is the proximity of the estimated parameters to the real ones and the sensitivity which tells us if we are able to distinguish between two different tumors with close parameters. The second analysis focuses on the shape of the minimization surface around the estimated point. The parameter estimation method, as explained in the previous section, minimizes the objective function C. The shape of this function around its minimum shows us the feasibility of the minimization process.

1) Analyses of the Estimated Parameters: In Figures 5(a) we show the estimated diffusion parameters along with the real ones. In order to demonstrate the results, we project the high dimensional parameter space onto the 2D (d w , d g ). The larger markers in the plot represent the real parameters used to grow the synthetic tumors and the smaller ones represent the estimated parameters retrieved from the images. Each small marker with a specific shape and color is the estimate for the larger marker with the same shape and color. Although there is only one estimate for each parameter set d w , d g , ρ there are multiple small markers for each large marker due to projecting onto lower dimensional space. In other words, different small markers of the same shape and color are the estimated parameters of the tumors with different ρ but same d w and d g .

Analyzing Figure 5 we observe that the parameter estimation method is able to retrieve the value of d w with good accuracy. Moreover, the method is able to distinguish between different tumors with close diffusion coefficients. The estimation of d g on the other hand seems to be less accurate. We notice the consistent positive bias in the estimate of d g which increases with increasing d w . We believe there are two reasons for this. The first one is the difference between numerical schemes we use to solve the reaction-diffusion PDE and the traveling time formulation. The numerical scheme for the PDE [START_REF] Mccorquodale | A cartesian grid embedded boundary method for the heat equation in irregular domains[END_REF] uses linear interpolation of the diffusion tensors between voxels creating higher diffusion within the grey matter neighboring white matter. The traveling time formulation, which uses the diffusion tensors on the voxels, accounts for this by increasing d g therefore estimating a higher d g . As a result as the value of d w increases the bias on d g increases. The second reason is computing the curvature effect term in Equation 9 using the images, where the contour enclosing the tumor delineation has sharp corners (due to discretization) which causes high curvature. Since the high curvature slows down the evolution, the traveling time formulation accounts for this by increasing the diffusion coefficient. This second reason is especially observed for the tumors where d w is low. Even in the presence of this bias we notice that for slowly diffusing tumors the d g estimates are very close to the real values and the method is able to distinguish between different tumors with close diffusion coefficients. For highly diffusing tumors the d g estimates are rather unreliable however the order of the ratio between d w and d g is well captured. The synthetic tumors are grown with the reaction-diffusion model with known parameters and synthetic images were created from these tumors. The parameter estimation method was applied to these images to retrieve the parameters of the model. The plot (a) shows the real diffusion rates dw and dg (the large markers) and the estimated diffusion rates (the small markers). Small markers of a specific shape and color are the estimates of the larger marker with the same shape and the same color. Regarding the estimation of T 0 , in Figure 5(b) we plot the estimated value of T 0 in the y-axis versus its real value in the x-axis where the y = x line is also drawn. Observing this plot we notice that the estimates for T 0 remains within the 10-15% margin of the real value, which shows that the proposed method is able to retrieve T 0 .

2) Analyses of the Minimization Surface: Regarding the shape of the minimization surface on the global scale, in our experiments we observed that this surface remains convex for all the tumors. However, the exact shape of the surface and its slope in different directions around the minimum point varied. We know that the estimated parameters provide us the best fit to the evolution of the tumor delineation observed in a set of images, let us say with an error of C * . The question we want to answer is how much this evolution varies from the optimum when we slightly move away from the "best" parameter set.

In order to answer this question, for an estimated parameter set (d * w , d * g , T * 0 ) that gives the minimum error of C * we find the other parameter sets which give an error smaller than C * + ǫ. In other words, parameter sets which provides an evolution of the tumor delineation which is ǫ away from the best fit in the average. In our high dimensional parameter space these parameter sets are enclosed in an ellipsoid around the estimated point which we name ǫ-ellipsoid (see Appendix D for details on how we construct the ǫ-ellipsoids). In Figures 6(a,b) we show the projections of some of these ǫ-ellipsoids on the respective parameter spaces where the round dots are the actual parameters, the crosses are the estimated parameters and ellipses around each cross are the projections of the ǫ-ellipsoids. In this study we have chosen to set ǫ = 0.2 mm 2 . Observing Figure 6(a) we notice that the major axis of the ellipses remain parallel to d g axis however, this is due to the difference of scale between d axis and the d g axis. When placed on the same scale these ellipses are rather circular. The second thing we notice is that the ellipses grow with increasing d w . This is a consequence of using normalized distances between surfaces in our error measure, see Equation 12. As d w increases the tumor diffuses faster in the white matter and its size increases. As a result the boundaries of the visible tumor reaches the extent of the white matter and most of the surface enclosing the tumor delineation in the image remains in the grey matter (as grey matter diffusion is much lower the tumor stops in the white-grey matter boundary) or reaches the boundaries of the brain. Therefore changing d w does not affect these portions of the surface and its contribution to the error measure decreases resulting in the larger ellipses we observe. This shows us that for more diffusive tumors a larger set of parameters yields similar errors therefore minimization surface is flatter. In Figure 6(b) we observe the coupling between d w and T 0 . One can obtain a similar evolution by increasing d w and decreasing T 0 (and vice-versa). The reason for this can be explained by the effect of convergence given in Equation 5, see Figure 2. We see that when T 0 is lower the speed of the tumor delineation is slower but if we increase the value of the diffusion we would obtain a similar evolution. The shape of the convergence curve in Figure 2 allows us to distinguish between these different cases and therefore find a minimum. In Equation 5we also notice that if T 0 is very high then a small change in T 0 does not affect the speed of the tumor delineation and this is the reason why we observe ellipsoids with major-axis parallel to the T 0 axis at high T 0 values. One can think of the extreme case where T 0 is very large and the effect of convergence becomes negligible. In this case we would expect its value not to change anything however, including the size of the tumor in the first image using the error term C 2 (Equation 16) helps us distinguish between very high T 0 values. One important conclusion we can reach from the sensitivity analysis presented in this section is the dependence of the estimated parameters on the uncertainty on the extracted tumor delineations. The tumor delineation either done by the expert manually or done using a segmentation algorithm has a variability. Kaus et al. in [START_REF] Kaus | Automated segmentation of MR images of brain tumors[END_REF] studied this variation for manual delineations and found out that in the case of low-grade gliomas the intra-expert variability is around 2% and the inter-expert variability is around 10% by volume comparison (they compared the volume of the tumor delineated by different experts and by the same expert multiple times). These variability values were greater than the automatic algorithms [START_REF] Kaus | Automated segmentation of MR images of brain tumors[END_REF]. The ǫ = 0.2 mm 2 value we used in this section corresponds to the inter-expert variation for manual delineation, 10% of volume variation for a tumor of 2.0 cm in diameter, see Appendix E for details. In other words, the ǫ-ellipsoids drawn in Figures 6(a) and (b) also demonstrates the influence of the inter-expert variability in the tumor delineations on the estimated parameters.

C. Varying the fixed ρ and Speed of Growth

In all the above experiments we have fixed the value of ρ to its real value. Naturally the diffusion rate estimates depend on this value of ρ. By fixing ρ we actually determine the location of the d w and d g estimates. In order to understand the effect of the value of ρ on the estimation of diffusion rates and the coupling between ρ and D, we have performed a slightly different experiment. Instead of fixing ρ to its real value we have set it to a different value and then estimated the other parameters d w , d g and T 0 . For the ease of demonstration we only show the estimation results for the synthetic tumors with ρ = 0.012/day. The experiment we performed is the same as the one explained in the previous section however, this time in q d (w,g) ρ values estimated by fixing ρ = 0.015 versus ρ = 0.012. We also plot the y = x line for a better comparison. We know that the estimated dw and dg values depend on where we fix the ρ. However, observing these figures we note that no matter what value we fix ρ to, the product of ρ and the estimated diffusion coefficient d (w,g) remains constant. Therefore the asymptotic speed of growth of the tumor in the white matter and in the grey matter can be estimated uniquely regardless of the choice of ρ.

the estimation method we set ρ = 0.015/day. As expected the estimated diffusion rates are lower than the values estimated by setting ρ = 0.012/day. The interesting point however, was not the change in the values but the coupling between D and ρ. In Figure 7(a) we plot v w = 2 √ d w ρ computed with ρ = 0.015 and the d w value estimated by fixing ρ to this value versus v computed using ρ = 0.012 and the d w estimated with this ρ. Figure 7(b) is the same plot for d g values. We observe from these graphs that the estimated diffusion rates change when we change the fixed ρ however, the product of the proliferation and the diffusion rates remain constant. The value v w,g = 2 d (w,g) ρ is the asymptotic speed of tumor growth and the proposed method is able to retrieve this speed uniquely from time series of images for all the 180 synthetic cases used in this analysis.

Here we would like to draw the attention to one important conclusion that can be reached with the observations and the experiments presented so far. In Figure 4 and in the accompanying table we have demonstrated that two tumors showing very similar volume changes in time had very similar asymptotic speeds in the white and in the grey matter (2 ρd w,g ). In this section, we have shown that using the parameter estimation method we can uniquely identify this product for all the tumors presented in this study. These two suggests that the evolution of the tumor delineation -for the synthetic tumors grown by reaction-diffusion models -can be uniquely matched by the proposed parameter estimation method and reproduced by the traveling time formulation. In Figure 8 we demonstrate this for a synthetic tumor. We first estimate the diffusion rates of the synthetic tumor whose evolution is shown in Figure 8(a). The estimation process is done once by setting ρ = 0.012/day and another time by setting ρ = 0.015/day. We have used three images in the estimation each 200 days apart. The resulting diffusion parameters are given in the table accompanying the images and the resulting optimum evolutions are shown in Figure 8(b) in red and green contours respectively. Following this, we start from the last image of the tumor used in the estimation and predict the further evolution of the tumor using the estimated parameters and the traveling time formulation. We do this once using the parameters estimated by setting ρ = 0.012 and once for those estimated by setting ρ = 0.015. We compare the prediction results with the actual evolution of the synthetic tumor in Figures 8(c) and (d) respectively. We observe that the predictions obtained by using different ρ values are almost identical and they show very high resemblance with the actual evolution of the tumor.

After analyzing the presented methodology with synthetic cases we apply the method to a few real cases in the next section.

D. Case Studies with Patient Data

The evaluation of parameter estimation for tumor growth models using real patient images is not easy because we do not have access to the real values of the parameters. The real values could be found using microscopic in-vivo analysis however, up to the best of our knowledge such a study has not been performed yet. In this work we perform an indirect evaluation for the proposed parameter estimation method using patient images. For a given patient dataset, we estimate the parameters using all but the image taken at the last time point. Then using the estimated parameters we simulate the evolution of the tumor delineation starting from the image taken just before the last one for the same number of days as the time difference between the last image and the one before it. We then compare the evolution predicted using the estimated parameters and the traveling time formulation with the one observed in the last image. The correlation between the prediction and the observed delineation provides us with a qualitative evaluation of the estimated parameters. The strong assumption we do here is that the values of the parameters remain constant between the images. Considering therapy and other effects on the tumor this assumption is not very realistic. However, we consider the estimated parameters as the average parameters over time including all the effects and carry on with the analysis.

As a preliminary step in this work we use two patient datasets which include anatomical and diffusion tensor MR images. The dataset for the first patient, who suffers from a high grade glioma (Glioblastoma Multiforme), includes T1-post gadolinium MR images (with the resolution of 0.5x0.5x6.5 mm 3 ) at three successive different time points and diffusion tensor MR image (with the resolution of 2.5x2.5x2.5 mm 3 ) taken at the second time point. The second patient suffers from a low grade glioma (second grade astrocytoma) and the dataset for this patient includes T2 flair MR images (with the resolution of 0.5x0.5x6.5 mm 3 ) at 5 successive time points and the DT-MRI image (with the resolution of 2.5x2.5x2.5 mm 3 ) taken at the first time point. For both cases the tumor boundaries were manually delineated by an expert in each image separately. We note that although manual delineations were used, in terms of the method and the analysis any segmentation algorithm can replace the manual delineation. As explained in the previous paragraph we estimate the parameters of the tumor growth model using all the images but the last one. The DT-MRI images of the patients are used to construct the diffusion tensor D of the tumor growth model. In constructing the diffusion tensor D for tumor cells we adapt the proposed models in [START_REF] Clatz | Realistic simulation of the 3d growth of brain tumors in mr images coupling diffusion with biomechanical deformation[END_REF] and [START_REF] Jbabdi | Simulation of anisotropic growth of low-grade gliomas using diffusion tensor imaging[END_REF]. As our focus in this article is the parameter estimation method we do not introduce a new diffusion tensor construction and we use the existing ones. Clatz et al. have proposed the tensor construction as given in Equation 3 for the high grade gliomas, following this we use this type of construction for our high grade case. On the other hand, Jbabdi et al. proposed to use another construction for the low grade gliomas given as

D(x) = d g I V (x) [diag(αe 1 (x)d w , d g , d g )] V (x) T , (18)
where V (x) is the eigenvector matrix obtained by decomposing the water diffusion tensor D water , e 1 (x) is the principal eigenvalue of the same tensor and α here is a normalization factor such that highest e 1 value in the brain becomes 1. The difference between this construction and the one given in Equation 3is that in this one tumor cells are assumed to diffuse much faster along the fiber and they diffuse very slowly in the transverse direction. In the construction the diffusion rate in the grey matter is used also for this transverse diffusion rate. As a result of such a construction the evolution obtained is much more anisotropic and creates more "spiky" tumors. Following the model assumptions made by the authors in [START_REF] Jbabdi | Simulation of anisotropic growth of low-grade gliomas using diffusion tensor imaging[END_REF] we use this type of construction for the low grade case. The images used to estimate parameters, the estimated parameters and the predicted evolution of the tumor delineations along with the real delineations are given in Figures 9 and10. In the images in both Figures, first we show the anatomical images at the time of detection and the intermediate images used in the parameter estimation. On the intermediate images we also plot the manual delineations for the underlying image (white contour) and the simulated evolution of the tumor delineation with the estimated parameters (dark contour) obtained in the course of estimation. Following this we start from the last image (in time) used in the parameter estimation and predict the evolution of the tumor delineation until the acquisition of the final image (which was not used in the estimation). In the corresponding images we show the anatomical MR image taken at the last time point showing the final state of the tumor along the tumor delineation predicted using the estimated parameters drawn as the dark contour. In the accompanying tables we provide the values of the estimated parameters.

In the images of the first patient in Figure 9 (a,b,c) the tumor showed evolution in two different regions. In the first region seen on the upper left corner of the images the tumor has a much larger volume, contains a necrotic core and exerts a visible mass effect. The second region on the other hand is a newly emerging lesion with no observable mass effect. This part is believed to be a diffused branch of the larger region however, no connection was visible in the images most probably due to slice spacing. We apply our analysis to the newly emerging part because it does not exert a mass effect and it is ideal for our analysis. Following the discussions given in Section III-A we fix the value of ρ to be able to estimate the diffusion parameters. The proliferation rate was set at ρ = 0.05/day around the suggested value in the literature [START_REF] Swanson | Virtual brain tumours (gliomas) enhance the reality of medical imaging and highlight inadequacies of current therapy[END_REF]. We observe in the image (c) that the prediction of the tumor delineation is in very good agreement with the final state of the tumor. This shows us that although we can provide the speed of growth only, the parameters combined with the growth model are in good agreement with the evolution of the tumor. Moreover, for the high grade tumor we estimated the speed of growth along the white matter as v w = 0.31 mm/day and in the grey matter as v g = 0.02 mm/day which are in good agreement with the literature [START_REF] Swanson | Virtual brain tumours (gliomas) enhance the reality of medical imaging and highlight inadequacies of current therapy[END_REF], [START_REF] Clatz | Realistic simulation of the 3d 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].

In the case of the low grade tumor shown in Figure 10, the correlation between the predicted tumor delineation and the final state of the tumor is in line with our previous arguments. We observe that the slow evolution of the tumor is well captured by the estimated parameters. For the proliferation rate we pick a lower value than the one in the previous case since it is a lower grade tumor. It was set to ρ = 0.008/day. Through the estimated diffusion rates we find the speed of growth along the white matter as v w = 0.08 mm/day and in the grey matter as v g = 0.004 mm/day.

Comparing the speed values estimated for the high grade and the low grade glioma we observe the expected difference. However, we would like to note that this difference is affected by the difference in the tensor construction method. In our experiments we have observed that when the diffusion rates for the same patient are estimated using the tensor construction given in Equation 18, the resulting values are higher. This is consistent with the fact that the tensors constructed with this method have lower diffusion in the transversal direction of the fibers. As a result a higher diffusion rate is needed to explain the same amount of growth. If we would like to compare the speed values for these two tumors we should keep this effect in mind. Therefore, the difference between speed of growth of these tumors are higher than the difference given between the above mentioned values.

IV. CONCLUSIONS

In this work we proposed and analyzed a parameter estimation method for the reaction-diffusion tumor growth models in the context of brain gliomas. The proposed methodology uses the evolution of the tumor, visible in the series of patient images, for estimating the parameters of the growth model. Using the patient images, the real 3D geometry of the brain and the tumor, tissue inhomogeneities and different diffusion properties are taken into account by the method. Moreover, unlike the previous methods that use the tumor cell density distribution, which is not available in the images, the proposed method formulates the evolution of the tumor delineation based on the reaction-diffusion dynamics. Such a formulation provides us a consistent framework in which the observables are the same as the model outputs and this removes the need of assuming a tumor cell density distribution in the images. To the best of our knowledge, this article constitutes the first work focusing on the automatic parameter estimation problem for reaction-diffusion tumor growth models using series of medical images and the real geometry of the patient.

In order to understand the theoretical properties of the underlying parameter estimation problem and the proposed method we performed thorough analysis using synthetic tumors for which the growth model parameters are known. The reaction-diffusion model that we have focused on has 3 different parameters: the diffusion rate in the grey matter d g , the diffusion rate in the white matter d w and the proliferation rate of tumor cells ρ. In our analysis we have shown that these parameters are coupled and therefore there is not a unique solution constrained by the observations made on medical images. However, we have also shown that once the proliferation rate ρ is fixed, then we were able to uniquely estimate the diffusion rates in the grey matter d g and in the white matter d w for all the 180 synthetic tumors presented in this work. Moreover, in this case we could also estimate the time elapsed between the emergence of the tumor and its detection, T 0 . In fixing ρ we assumed that its value can be found through microscopic analysis of biopsy results [START_REF] Schröder | The relationship between Ki-67 labeling and mitotic index in gliomas and meningiomas: demonstration of the variability of the intermitotic cycle time[END_REF], [START_REF] Johannessen | The clinical value of Ki-67/MIB-1 labeling index in human astrocytomas[END_REF]. We have also shown that the value of ρ determines the estimates of the other parameters. Investigating the coupling between diffusion and the proliferation rate we have shown that no matter what ρ value we fix, the product of the estimated diffusion rates with ρ remains constant for the same tumor. This product represents the speed of growth of the tumor in terms of the model parameters. We have shown that for each tumor used in this work we were able to uniquely estimate this product in the white matter and in the grey matter separately. Moreover, we have also shown that the estimated diffusion rates and the fixed proliferation rate were able to capture the growth dynamics of the tumor and simulate its evolution regardless of the value ρ was fixed to. This demonstrated that the proposed methodology is successful in creating the "patient-specific" model and perform personalized simulations.

We also applied our method to two real cases, one high grade glioma and one low grade. In the light of the synthetic analysis, we have set the proliferation rates ρ for these tumors to average values and estimated the diffusion rates. Using the estimated parameters and the traveling time formulation we have shown promising preliminary results in personalizing reaction-diffusion models. The prediction studies provided us an indirect validation for the estimated parameters and the formulation. Moreover, these results demonstrated that although the complete parameter estimation problem has a non-unique solution, by fixing ρ and estimating the diffusion rates we are able to personalize the growth model and use it to simulate the tumor growth and predict the further evolution. Therefore, the proposed methodology showed itself to be a successful attempt for adapting the tumor growth models to patient images and creating "patient-specific" models. The strongest assumption we made during the analysis of the real cases was that the parameters of the growth model do not change in time and they do not vary in space. This may not be realistic for the exact values of the parameters considering the existence of different types of therapies and the random nature of the tumor progression, which exists both in time and in space. However, we have regarded the estimated parameters as the average values over time and space including all the different effects, which is a clinically logical step [START_REF] Mandonnet | Continuous growth of mean tumor diameter in a subset of grade II gliomas[END_REF]. On the other hand, independent parameter estimation and analysis could be done between each set of two successive images as well. Such an analysis combined with the time course of the therapy could give us hints on the effect of the therapy on different parameters and on the growth speed of the tumor.

The method proposed in this article is a first attempt to solve the parameter estimation problem and there are different improvements that should be integrated in the future studies. As a first step, we ignored the mass effect of the tumor. In most glioma cases the mass effect is apparent, smaller in the low grade gliomas and larger for the higher grades. For a complete modeling, the mass effect should be taken into account in the parameter estimation methodology. The second point is the effect of the tumor growth on the fiber structure. The white matter is disturbed due to the tumor invasion and this should be included to have a more complete methodology [START_REF] Price | Diffusion tensor imaging of brain tumours at 3t: A potential tool for assessing white matter tract invasion?[END_REF]. The third point for further improvement is the therapy response and the shrinkage of the tumor. In this work we focused on the reaction-diffusion growth models which only formulate the growth of the tumor. Therefore, the shrinkage due to therapy or any other process is not taken into account in our methodology. Fortunately, the proposed traveling time formulation and the parameter estimation methodology can be used in the case of shrinkage after certain modifications (embedding the Eikonal formulation in a level-set framework). However, the parameter estimation problem should be redefined to decouple the therapy response and natural growth of the tumor. In the further studies the effect of the therapy and the possible shrinkage of the tumor should be taken into account. One other point that can be included in the methodology is the multi-phase modeling. Multi-phase models describe and formulate the interactions between different phases such as different tumor populations and tumor-brain interface. We can think of generalizing the proposed method to include multi-phases through coupling several Eikonal models. Such an approach would enable us to easily explain multi-foci gliomas and also model the response of the brain tissue to tumor growth. As a last point for improvement, we aimed not to use any assumptions on the tumor cell density in our method however, as more information becomes available on the cell density distribution, it can also be integrated into the proposed methodology.

In this work our aim was to propose the parameter estimation method and analyse it theoretically. We also applied the method to a few real cases demonstrating the applicability of the method and showed promising results. Eventually a more thorough analysis of the estimated parameters and the estimation methodology should be performed using a big dataset of patient images. Such a study will let us better understand the clinical signficance of the estimated parameters and the constructed patient-specific models. In the follow-up of this work we plan to focus on this direction. There are several problems that should be overcome for this purpose. The first problem is the lack of diffusion tensor imaging for the patients. As we have seen the DTI is very important in the modeling and in the estimation of the parameters therefore, it is crucial to have this information. The advances in the registration methods can be helpful to solve this problem as they would give us the opportunity to create and register DT-MRI atlases on the patient images. The second problem is regarding the surgery applied in glioma cases. The surgery changes the structure of brain as well as the properties of the tumor. In order to overcome this problem, we need to adjust the traveling time formulation such that it can describe the evolution of the tumor delineation between pre-op and post-op images.

The method proposed in this article adapts the reaction-diffusion model to specific patient images. There are two clinically relevant outputs of this adaptation: the estimated parameters and the "patient-specific" growth model. The estimated parameters, specifically the speed of growth and the measure of differential motility, can serve as quantification measures for tumor growth and help the diagnosis process. Through population studies and group analysis, correspondence between different parameter ranges and different WHO grades of gliomas can be found [START_REF] Mandonnet | Continuous growth of mean tumor diameter in a subset of grade II gliomas[END_REF]. This analysis may provide us with information about variation within tumor grades and could allow to identify the transition from low grade glioma -which is not treated but only monitored -to a high grade tumor -requiring immediate treatment -more accurately. A quantitative analysis of the tumor evolution during therapy, i.e. an accurate, localized estimation of tumor growth or shrinkage, may allow to estimate the efficacy of a chosen treatment option much quicker than with current qualitative approaches. Determining the changes in growth parameters due to the treatment may finally provide the means to understand the tumor response to therapy -specific for the individual patient [START_REF] Swanson | Dynamics of a model for brain tumors reveals a small window for therapeutic intervention[END_REF]. The patient-specific tumor growth model, in other words the generic model with the patient-specific parameters, gives us the opportunity to simulate the specific evolution of a patients tumor. The expected anisotropic growth suggested by the patient-specific model may guide biopsies or at least more targeted imaging methods (such as MRSI) for more accurate estimations of the tumor extensions. Eventually, as the generic models become more realistic patient-specific models can be used to better plan the therapy process [START_REF] Konukoglu | Extrapolating tumor invasion margins for physiologically determined radiotherapy regions[END_REF] and predict possible outcomes.

APPENDIX A THE TRAVELING TIME FORMULATION

The equation of the approximation for the speed of the tumor front including the effect of time convergence is given as

v(t) = √ n ′ Dn(2 √ ρ -3/(2t √ ρ)). (19)
This equation is valid when the front of the tumor is not curved and the parameters are constant. In more general case, as in the case of the brain, the front is curved and the parameters vary. In order to formulate the motion of the tumor front we make a voxel based assumption. We assume that within the voxel the tumor front is planar and the parameters of the model (D and ρ) are constant and the values of them are taken as the values at that voxel. Under this assumption Equation 19 can be converted into a traveling time formulation for the tumor front using the same idea as explained in [START_REF] Sethian | Level set methods and fast marching methods: Evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science[END_REF]. Using the relation

|∇T | = 1/v(t) = 2 ρn t Dn - 3 2T n t Dn ρ -1 (20)
where T is an implicit time function such that it embeds the location of the tumor delineation at different times as iso-time surfaces. Remembering that the n can be written as ∇T /|∇T | (because T is an implicit function) we can write the traveling time formulation as

√ ∇T ′ D∇T = 2 √ ρT 4ρT -3 . (21)
This equation alone only gives the relation of successive iso-time surfaces of T . In order to build the solution we need a Dirichlet type boundary conditions, namely an initial surface for which we know the T value. In the context of the tumor growth modeling this surface is given as the surface enclosing the tumor delineation in the image. Using the image we can write the necessary Dirichlet condition as

T (x) = T 0 ∀x ∈ Γ (22)
where T 0 is the initial time we start from and Γ is the surface found in the image.

APPENDIX B DERIVATION OF THE EFFECT OF THE CURVATURE FOR ANISOTROPIC DIFFUSION TENSORS

This derivation follows the derivation given in [START_REF] Keener | Mathematical physiology[END_REF] and modifies it for the anisotropic tensor case. The reaction-diffusion model has the general form:

u t = ∇ • (D∇u) + ρu(1 -u). (23)
We apply a coordinate change by parameterizing the moving frame of the u function as

x = X(ξ, τ ), t = τ. (24)
We assume that this parameterization is a diffeomorphism. By chain rule the partial derivatives using the new coordinates can be written as

∂ ∂ξ i = ∂X j ∂ξ i ∂ ∂x j (25)
∂ ∂τ = ∂ ∂t + ∂X j ∂τ ∂ ∂x j (26)
where the indices are summed (this is the case throughout this appendix). Likewise the partial derivatives with respect to the Euclidean coordinates can be written in terms of the new coordinate system.

∂ ∂x i = α ij ∂ ∂ξ j (27)
∂ ∂t = ∂ ∂τ -α jk ∂X j ∂t ∂ ∂ξ k (28)
represents are the partial derivatives in terms of the new coordinate system. α ij is the ij th component of the inverse of the Jacobian matrix with respect to the parameterization X. We identify ξ 1 as the normal direction to the iso-surfaces of u at every point. We also define the tangent and the normal vectors of the parameterization as

r i = ∂X j ∂ξ i (29)
n i = r j × r k , j, k = i. (30)
Using this we can define the [α] matrix using these vectors:

α ij = (n j ) i r j n j . (31)
For the ease of derivation, through the choice of the parameterization we let r 1 • r 2 = 0 and r 1 • r 3 = 0 (r 1 n 1 ) and set r 1 as the normal vector to the iso-surface of the u function. The derivative terms in the reaction-diffusion equation become

u t = ∂u ∂τ -α jk ∂X j ∂τ ∂u ∂ξ k (32) ∇ • (D∇u) = α kp α ij d ki ∂ 2 u ∂ξ p ∂ξ j (33)
+ ∂ ∂x k (d ki α ij ) ∂u ∂ξ j .
Then the whole equation can be written as

α kp α ij d ki ∂ 2 u ∂ξ p ∂ξ j + ∂ ∂x k (d ki α ij ) ∂u ∂ξ j (34)
-( ∂u ∂τ -α jk ∂X j ∂τ ∂u ∂ξ k ) + ρu(1 -u) = 0.
At this moment we use the two big assumptions made in [START_REF] Keener | Mathematical physiology[END_REF]. The first of these assumptions say that the spatial variation of ξ 1 is much smaller than the ξ 2 and ξ 3 . This actually says that the normal to the front changes faster than the tangent space of the parameterization. Therefore, the effect of curvature is in a lower order than the speed of the moving frame. Remembering that the [α] is the inverse Jacobian matrix of the parameterization X this assumption lets us say that α j1 = O(1) while α jk = O(ǫ). The second big assumption is that to the leading order in ǫ, u is independent of τ . In the planar evolution this assumption readily holds since the solution of the reaction-diffusion equation is a traveling wave therefore, does not depend on time. However, for the curved evolution this does not have hold. This assumption on the dependence on τ lets us treat the curved evolution as if it admits a traveling wave. Using the singular perturbation method we can gather the first order terms and 34 reduces to

α k1 d k1 α i1 ∂ 2 u ∂ξ 2 1 + ∂ ∂x k (d ki α i1 ) ∂u ∂ξ 1 (35)
α j1 ∂X j ∂τ ∂u ∂ξ 1 + ρu(1 -u) = O(ǫ).
Gathering the terms and recognizing the matrix multiplications this equation can be rewritten in the compact form

α ′ Dα ∂ 2 u ∂ξ 2 + (∇ • (Dα) + α • ∂X ∂τ ) ∂u ∂ξ 1 + ρu(1 -u) = O(ǫ). (36)
Now in order to have traveling wave solution this ODE should have the same form as the one in Equation 4. This means that we need the coefficients of this equation to be constants. However, this will not be possible for every iso-contour of the function u. The curvature will have different effects for different iso-surfaces. Hence, we require it only for the origin of the moving frame (u = 0.5 iso-contour in the case of logistic growth). Using this information and the restriction we obtain

α ′ Dα = ρ (37) ∇ • (Dα) + α • ∂X ∂t = 2ρ. (38)
At this point we remember that α n 1 which is normal to the iso-surface of u. We define a level set function S such that the zero-level set of S will correspond to the origin of our moving frame therefore, ∇S/|∇S| = n. We can then write α as α = -K∇S where K is just a coefficient to be determined. From Equation 37 we find K as

K = ρ ∇S ′ D∇S . (39)
On the other hand, the Equation 38gives us

-∇ • (DK∇S) -K∇S • X t = ρc 0 . (40)
In order to replace X t we need one more relation which comes from the fact that the value of function S on the origin of the moving frame doesn't change by construction. Therefore, ∂ ∂t S(x, t)| on the moving frame origin = 0 (41)

∇S • X t + S t = 0. (42)
Placing this in Equation 40we obtain

∇ • (DK∇S) + KS t = ρc 0 (43) ∇ • (D∇S • ρ ∇S ′ D∇S ) + ρ ∇S ′ D∇S S t = ρc 0 . (44)
Now transforming the dynamic Hamilton-Jacobi equation given above into a static one by inverting the embedding method explained in [START_REF] Osher | A level set formulation for the solution of the dirichlet problem for hamilton-jacobi equations[END_REF] we obtain the anisotropic Eikonal equation with the curvature term

∇ • ((D∇T ) ρ ∇T ′ D∇T ) + ρ ∇T ′ D∇T = ρc 0 . (45)
Relocating terms and letting c 0 = 2 we get our formulation:

[2 √ ρ -∇ • ( D∇T √ ∇T ′ D∇T )] √ ∇T ′ D∇T = 1. (46)
This equation uses the asymptotic speed v given in Section II-A. We can also replace v by v(t) and obtain the whole equation

{ 4ρT -3 2 √ ρT -∇ • D∇T √ ∇T ′ D∇T } √ ∇T ′ D∇T = 1 (47)

APPENDIX C NUMERICAL METHOD FOR THE ANISOTROPIC EIKONAL EQUATION WITH CURVATURE EFFECT

The static Hamilton Jacobi equation given in Equation 5 is a first order equation and has the form of an anisotropic Eikonal equation:

F (x) √ ∇T ′ D∇T = 1, (48)
where the additional F (x) is a spatially varying speed function. In [START_REF] Konukoglu | A recursive anisotropic fast marching approach to reaction diffusion equation: Application to tumor growth modeling[END_REF] we have proposed a numerical method to solve this kind of equation. The proposed method is based on the Fast Marching methods [START_REF] Sethian | Level set methods and fast marching methods: Evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science[END_REF] and modifies it in order to take into account the anisotropy in the equation. The original Fast Marching method numerically solves the isotropic Eikonal equation (F (x)|∇T | = 1) by following the gradient direction of T as it constructs the solution. This, when applied to the Equation 48, creates erroneous solutions because the characteristic direction should be followed when constructing the solution of this anisotropic equation [START_REF] Sethian | Ordered upwind methods for static hamilton-jacobi equations: theory and algorithms[END_REF].

The numerical method proposed in [START_REF] Konukoglu | A recursive anisotropic fast marching approach to reaction diffusion equation: Application to tumor growth modeling[END_REF] integrates a recursive correction scheme inside the original fast marching algorithm. This correction scheme makes sure that as we solve the Equation 5 we follow the characteristic directions of the PDE [START_REF] Kevorkian | Partial differential equations: Analytical solution techniques[END_REF]. The overall method starts from a given initial contour and sweeps the domain outwards finding the solution of the Equation 48at each voxel. This equation has 2 different solutions at each voxel and in the case of Equation 5we choose the value such that as we move away from the delineation the T value increases (since the tumor delineation will pass from those points at a later time in the case where the tumor grows). The advantages of this numerical method are that it is a sweeping method and it only uses the immediate neighbors of a point rather than using points far away [START_REF] Sethian | Ordered upwind methods for static hamilton-jacobi equations: theory and algorithms[END_REF] to compute the values. Therefore, it is a fast and accurate method for solving anisotropic Eikonal equations. For the details of the algorithm please refer to [START_REF] Konukoglu | A recursive anisotropic fast marching approach to reaction diffusion equation: Application to tumor growth modeling[END_REF]. Here, regardless of the details of the algorithm, we continue our discussion based on the fact that we have a sweeping algorithm which solves anisotropic Eikonal equations in a fast manner.

Using the sweeping method for solving the Equation 7(and for the Equation 8) is not very obvious because it is not a first order equation due to the divergence term. These equations can be solved with other iterative methods [START_REF] Osher | A level set formulation for the solution of the dirichlet problem for hamilton-jacobi equations[END_REF], [START_REF] Qian | A fast sweeping method for static convex hamilton-jacobi equations[END_REF] however, these methods are not very fast. In order to benefit from the advantages of the sweeping method explained above, we separate the curvature part from the equation and construct an iterative method that solves anisotropic Eikonal equations at each iterations with different speed terms. The form we use for Equation 7(it is the same construction for Equation 8) becomes

{ 4ρT -3 2 √ ρT -∇ • D∇T √ ∇T ′ D∇T } √ ∇T ′ D∇T = 1 (49)
{ 4ρT -3 2 √ ρT + F curv } √ ∇T ′ D∇T = 1. (50)
Viewing the convergence term as a speed term independent of T as F curv enables us to use the sweeping method and construct the simple iterative method

F 0 curv = 0 (51)
{ 4ρT n-1 -3 2 √ ρT n-1 + F n-1 curv } √ ∇T ′n-1 D∇T n-1 = 1 (52) Compute T n-1 (53)
C n-1 = -∇ • D∇T n-1 √ ∇T ′n-1 D∇T n-1 (54)

F n curv = F n-1 curv + α(C n-1 -F n-1 curv ). ( 55
)
where α n-1 < 1 is the parameter determining the rate of convergence which in our case is taken as α = 0.8. In Equation 55we see that the F n curv is updated with a proportional gain using the error made in the previous iteration. In this respect this scheme is similar to the feedback control loops. We iterate this algorithm until

x∈Ω |C n -F n curv | < ǫ (56)
where we sum represents the summation over all points in the domain of computation and ǫ is a small value. Once this criteria is satisfied we know that F curv is indeed the effect of the curvature. The rate of convergence depends on α however, in our experiments we have observed that for a large range of α ∈ (0.2, 0.8) the rate is pretty rapid see Figure 11.

D CONSTRUCTION OF ǫ-ELLIPSOIDS

The estimation problem in this work is formulated as the optimization problem with the objective function C. For a given set of images the method tries to find the parameters of the tumor growth model which would minimize the value of the function C. ǫ-ellipsoids is a simple way to understand the shape and the steepness of the minimization surface around the minimum point. The construction of the ǫ-ellipsoids is as follows. For a given parameter estimation problem let us say the estimated parameters p * corresponds to an error value of C * . As a consequence C * is the minimum of the objective function C for this problem. We first construct the quadratic approximation of C around p *

C = C * + g ′ (p -p * ) + 1 2 (p -p * ) ′ G(p -p * ), ( 57
)
where g is the gradient vector and G is the Hessian at p * . Since p * is the minimum we know that g = 0. Moreover since the point p * is the minimum of C the G is a positive definite matrix. The construction of the quadratic approximation is done by sampling the function C and fitting a quadratic function by least square minimization.

Once the quadratic approximation of C is obtained we define the ǫ-ellipsoid as follows

P = {p|C(p) = C * + ǫ}, (58)
where the set P is the ǫ-ellipsoid and p is an arbitrary parameter set. Since G is a positive definite matrix we are sure that P is a closed surface and for all the points remaining inside P , C(p) < C * + ǫ.

Using the ǫ-ellipsoid we enclose a set of parameter sets for which each parameter set produces an evolution of the tumor delineation that is ǫ close to the optimum evolution created by p * . This means if the ǫ-ellipsoid is big for a problem then the minimization surface is flatter therefore, it is harder to find the minimum point. Moreover, the directions of the semi-major and semi-minor axis of the ellipsoid provides us the coupling between different parameters.

E RELATIONSHIP BETWEEN ǫ-ELLIPSOIDS AND UNCERTAINTY ON TUMOR DELINEATIONS

The method presented in this work estimates the parameters of the reaction-diffusion tumor growth models using the tumor delineations extracted from the patient MR images. Although we do not focus on the segmentation algorithm in this work (and we note once again that any segmentation algorithm could be used with the proposed method), we are interested in the variability of tumor delineation extracted from the images. The tumor boundaries either extracted by an expert manually or by a segmentation algorithm automatically contains a certain variability. This variability has been studied by Kaus et al. in [48]. In this work, they show that in the case of low-grade gliomas the highest variability was seen between the delineations done by different experts. This inter-expert variability is shown to be around 10% by volume overlap for the dataset they used.

In Section III-B2 we presented the overall variation of the estimated parameters by showing the shape of the minimization surface around the optimum parameters. In doing so we have shown the extend of the basin around the optimum parameters using the ǫ-ellipsoids as explained in Appendix D. Here we show the relationship between the variability of input tumor delineations and the ǫ-ellipsoids. We show that the variability of the estimated parameters due to the uncertainty on the tumor delineations is given by the ǫ-ellipsoids. First we start by proving the following claim. Claim 1. Let Γ R1 and Γ R2 be two delineations of the same tumor in the same image obtained using different methods (same expert different time, different experts, different methods). Then the optimum parameters obtained using Γ R2 remains within the 2δ 2 -ellipsoid around the optimum parameters obtained using Γ R1 . And δ 2 = max (̺ (Γ R2 , Γ R1 ) , ̺ (Γ R1 , Γ R2 )), where ̺ is defined in Equation 14.

Proof: Let Γ S1 and Γ S2 be the simulated tumor delineations computed with the optimum parameteres obtained by using Γ R1 and Γ R2 respectively. Without loss of generality we assume that Γ R1 , Γ R2 , Γ S1 and Γ S2 are represented by the point sets {x n 1 }, {x n 2 }, {y n 1 } and {y n 2 } who have the same number of elements, M . Let us also define the functions

Fig. 1 .

 1 Fig. 1. MR Flair images of a grade II astrocytoma: (a) image at the first examination (b) image at the second examination. In the anatomical MR images we observe the evolution of the boundary of the visible part of the tumor rather than the tumor cell densities.

Fig. 2 .

 2 Fig.2. The reaction-diffusion equation in infinite cylinder admits a traveling wave solution resulting in a tumor front moving at constant speed. (a) The traveling wave has an asymptotic speed shown in dashed curve however, when we observe the speed of the iso-density contour u = 0.5 in time we notice the low rate of convergence to this speed (see the solid curve). An approximation of the speed of the iso-density contour including the convergence effect yields a closer curve to the the observed on (see point-dashed curve). (b) Starting from the same point the integrals of the speed curves, the distances to the initial point, are shown. Notice that we get a much better approximation when we add the convergence effect. All axis are in non-dimensional coordinates.

Fig. 3 .

 3 Fig. 3. Comparison between the reaction-diffusion model and the Eikonal Approximation:The temporal evolution of the iso-density contour is demonstrated for a synthetic tumor. Contours are shown for days 400, 600, 800, 1000 and 1200 from the innermost to outermost respectively. The synthetic tumor is virtually grown using the reaction-diffusion model. White contours are obtained by thresholding the tumor cell densities at u = 0.4 for the respective day values (400-600-800-1000-1200). Then in order to simulate the evolution of the iso-density contour (assumed to correspond to tumor delineation in real images) starting from day=400, without the knowledge of the tumor cell density distribution we use the traveling time formulation. Black curves are the contours we obtain at days 600 (2nd innermost) to 1200 (outermost). The average Euclidean distance between the black and the white contours is 0.78 ± 0.69 mm. We notice that the traveling time formulation is quite accurate in describing the evolution of the tumor delineation in the case of synthetic tumors. The tumors were grown in the images of a healthy subject for whom we also have the DT-MRIs. Parameters: (dw = 0.25 mm 2 /day,dg = 0.01 mm 2 /day,ρ = 0.012 day -1 )

Fig. 4 .

 4 Fig. 4. In the image we show the evolution of two different synthetic tumors virtually grown using the reaction-diffusion model with different parameters. The contours of the same color are the tumor delineations for the same tumor in 4 different images taken at 4 successive time instances (first image, 200, 300 and 400 days after the first image).The reaction-diffusion model parameters for these tumors are given in the table. We observe that although the diffusion and proliferation rates of these tumors are different the evolutions are almost the same. The difference between these evolutions measured using C is 0.644mm 2 which is lower than the minimum error we find by estimating the parameters using the traveling time formulation (C = 1.28 for red and C = 1.29 for green). This shows that we cannot distinguish between these two parameter sets if we observe either of the evolutions. We also observe that the products dwρ and dg ρ are very close for the two tumors. This tells us that although distinguishing between dw, dg and ρ is not obvious estimating the product of these values can be possible.

Fig. 5 .

 5 Fig.5. The results of the parameter estimation from time series of images for the synthetic tumor experiments. The synthetic tumors are grown with the reaction-diffusion model with known parameters and synthetic images were created from these tumors. The parameter estimation method was applied to these images to retrieve the parameters of the model. The plot (a) shows the real diffusion rates dw and dg (the large markers) and the estimated diffusion rates (the small markers). Small markers of a specific shape and color are the estimates of the larger marker with the same shape and the same color. Figure(b) plots the estimated initial time estimate T 0 (the time elapsed between the emergence of the tumor and the detection) vs. its real value. y = x line is also drawn for better comparison.

  Fig.5. The results of the parameter estimation from time series of images for the synthetic tumor experiments. The synthetic tumors are grown with the reaction-diffusion model with known parameters and synthetic images were created from these tumors. The parameter estimation method was applied to these images to retrieve the parameters of the model. The plot (a) shows the real diffusion rates dw and dg (the large markers) and the estimated diffusion rates (the small markers). Small markers of a specific shape and color are the estimates of the larger marker with the same shape and the same color. Figure(b) plots the estimated initial time estimate T 0 (the time elapsed between the emergence of the tumor and the detection) vs. its real value. y = x line is also drawn for better comparison.

Fig. 6 .

 6 Fig.6. Figures plot the projections of some of the ǫ-ellipsoids on the respective parameter spaces. The round dots are the real parameters of the reactiondiffusion model, the crosses are the estimated parameters and ellipses are the projections of the ǫ-ellipsoids for each cross. For a given cross, the cross represents the minimum of the respective minimization surface with an error of C * and all the points inside the ellipse surrounding that cross are the parameters who has error less than C * + 0.2 mm 2 . In other words ellipses enclose all the parameters producing very similar evolutions of the tumor delineation as the cross in the center.

vFig. 7 .

 7 Fig.[START_REF] Patel | A cellular automaton model of early tumor growth and invasion[END_REF]. In the figures we plot v (w,g) = 2 q d (w,g) ρ values estimated by fixing ρ = 0.015 versus ρ = 0.012. We also plot the y = x line for a better comparison. We know that the estimated dw and dg values depend on where we fix the ρ. However, observing these figures we note that no matter what value we fix ρ to, the product of ρ and the estimated diffusion coefficient d (w,g) remains constant. Therefore the asymptotic speed of growth of the tumor in the white matter and in the grey matter can be estimated uniquely regardless of the choice of ρ.

Fig. 8 .

 8 Fig. 8. The proposed parameter estimation methodology can uniquely estimate the asymptotic speeds of the tumor in the white and in the grey matter as a function of the parameters of the model. Figures show prediction experiments on synthetic tumors showing that combined with the traveling time formulation the estimated parameters capture the growth of the tumor regardless of the value the ρ is fixed to. (a) Evolution of the synthetic tumor shown in 3 white contours representing the delineation at the time of detection, 200 days and 400 days after the detection. (b) Optimum evolutions obtained by estimating the diffusion rates by setting ρ = 0.012/day in red and ρ = 0.015/day in green. The estimated parameters are given in the table. (c) Starting from the final image used in estimation (outermost white contour in (a)), further evolution is predicted using the parameters estimated when ρ was set to 0.012. The predicted evolution shown in red while the actual evolution shown in white. (d) Same image is shown for the prediction obtained using the parameters estimated when ρ was set to 0.015. Prediction shown in green contours. Observe that regardless of the fixed ρ value the traveling time formulation and the estimated parameters capture the growth of the tumor and able to simulate its evolution.

  1/day 0.66 mm 2 /day 0.0013 mm 2 /day

Fig. 9 .

 9 Fig. 9. The parameter estimation method is applied to the images of a real patient suffering from high grade glioma. Images in columns (a) and (b) shows different slices of the T1-post gadolinium images which are used to estimate the parameters of the growth model given in the table. In (b) we also show the manual delineation of the tumor (in white) used in parameter estimation along with the optimum simulation obtained by the estimated parameters (in black) (only white contour is shown in (a) since it is the same as the black one). (c) The final image showing the final state of the tumor and the evolution of the delineation predicted by the estimated parameters as the black contour.

  1/day 0.20 mm 2 /day 7x10 -4 mm 2 /day

Fig. 10 .

 10 Fig. 10. As a second case we applied our methodology to the images of a patient suffering from a low grade tumor. Images (a)-(d) show different slices of the T2 flair images and the manual delineations (in white) which are used to estimate the parameters of the growth model given in the table. Also in these images we show the simulated evolution of the tumor delineation obtained by the estimated parameters in black contours. The simulated evolution starts from the white contour in the Image (a). Images (e) are the slices of the final image showing the final state of the tumor and the delineation predicted by the estimated parameters as the black contour.

Fig. 11 .

 11 Fig. 11. The curve of convergence for the iterative scheme given by Equations 51-55. We visualize the difference P x∈Ω |C n -F n curv | as a function the iteration number n. The curve is obtained for the example shown in Figure 3.

This work has been partly supported by the European Health-e-Child project (IST-2004-027749), by the CompuTumor project (http://wwwsop.inria.fr/asclepios/projects/boston/) and by Microsoft Research, Cambridge.