Diffusive modelling of glioma evolution : a review

Gliomas, the most aggressive form of brain cancer, are known for their widespread invasion into the tissue near the tumor lesion. Exponential models, which have been widely used in other types of cancers, cannot be used for the simulation of tumor growth, due to the diffusive behavior of glioma. Diffusive models that have been proposed in the last two decades seem to better approximate the expansion of gliomas. This paper covers the history of glioma diffusive modelling, starting from the simplified initial model in 90s and describing how this have been enriched to take into account heterogenous brain tissue, anisotropic migration of glioma cells and adjustable proliferation rates. Especially, adjustable proliferation rates are very important for modelling therapy plans and personalising therapy to different patients.


INTRODUCTION
Glioma is a type of cancer of central nervous system that starts in the brain or in the spine.It is called a glioma because it arises from glial cells.The most common site of gliomas is the brain.Gliomas constitute more than 50% of all brain cancer cases.Despite extended research in this area, patients are rarely given more than 12 months survival time [1].
The diagnosis of gliomas can be done by MRI, CT, angiogram or biopsy, with MRI being the most common method.After glioma is diagnosed, treatment is directly necessary.However, the most important problem with diagnosis and treatment is that gliomas are characterized by invasiveness.This means that there are cells diffused beyond the imaged tumor, which cannot be visualized by common imaging techniques.Thus, even if the clinician allows some safety margin during resection around the imaged tumor, cancer is expected to recur.
These dispiriting results have forced researchers worldwide to work with understanding the glioma grow-th procedure and its special pathology.One of the rapidly emerging fields in glioma study is tumour growth modeling; researchers have been working on finding mathematical models that efficiently describe the glioma growth procedure.

Mathematical Modelling
A mathematical model uses mathematical language to describe a system.In the case of glioma, a mathematical description of how glioma grows is under research.At first someone has to bear in mind the basic features that any typical mathematical model in biology must possess [2].These are schematically described in Figure 1.Firstly, the model should be initiated within a realistic biological state.Additionally, the modeled biological processes should be understood and discretized as much as possibly, meaning that steps and real biological parameters should be isolated.Continuing, it is essential to allocate a mechanism that could simulate these steps and incorporates these parameters.Specifically, this mechanism could be described by an equation.Going further, the next step is to study the model mathematically and come up with solutions that include realistic boundary and initial conditions.Lastly, after having acquired the theoretical results it is of great importance to go back in biological process with predictions, comments and suggestions for experiments that will either ascertain or disprove the developed model.At this level, model success is highly dependent on combining experimentation and theory together.Because, even if the experimental results indicate that the model is incorrect, this is the right way to reach a successful conclusion.In final consideration, mathematics is very important in biology, however, they must be treated with seriousness.If mathematics is used for solving any biological process, without thoroughly studying the biological background, it is very possible to come up with solutions that not only do not contribute to corroborated conclusions, but also do harm.As stated in [2], the theoretical literature abounds with many such articles.
Recently, with the explosion of biological sequence data, many biological sequence databases have redun-dant sequences which can cause problems for data analysis.These redundant sequences cannot provide valuable information for analysis but detracts from the statistical significance of interesting hits.Moreover, processing these redundant sequences often requires more time and computational resources.Removing redundant sequences is undoubtedly very helpful for performing statistical analysis and accelerating extensive database searching [1].And it is also a way to obtain the real protein families and their representatives from a large sequences dataset.Therefore, it is necessary to develop an appropriate algorithm to remove redundant sequences from a biological sequence database.
Hobohm and Sander's algorithm is a widely used algorithm in many redundant sequence removing programs.Hobohm and Sander's algorithm was firstly introduced by U. Hobohm et al. of EMBL laboratory in 1992.In 1998, Lissa Holm and Chris Sander developed a program based on this algorithm to generate a nonredundant protein database NRDB90 [2].After that, other researchers developed some programs for removing redundant sequences on the basis of Hobohm and Sander's algorithm, such as CD-HIT and PISCES.
A lot of research is currently taking place in mathematically modeling the glioma growth procedure.An efficient mathematical model for gliomas could help researchers and clinicians to get a better understanding of tumour growth pathology.They could also to predict the aggressiveness of a patient's tumour and, thus, to better define the margins of a tumour, so as to use them when applying resection or radiotherapy.Moreover, by importing therapy parameters into the model, the clinician can predict which therapy scheme is expected to yield better results for the patient.

General Cancer Modelling
Studying tumour expansion and simulating this according to mathematical models, has been an area of studies in cancer since late 90s' [3][4][5].Tumour growth has been studied by a series of models.The first models focused on tumour behavior in time.More specifically, the first proposed temporal models were based on either exponential, logistic or Gompertz laws [6].As expected, these models were followed by spatial growth models in later years.Thus, one such deterministic model has been used to simulate cancer growth as a wave phenomenon, taking into account mitosis and nutritient depletion [7].Moreover, deterministic models taking into account immune response [8] or mitotic rates changes [9] have been proposed.
The research efforts for applying these models in infiltrative cancers, such as glioma, failed because cell motility hadn't been included in the model.Thus, some stochastic and cellular automata have been introduced for glioma simulation, taking into account cell cycle, lack or abundance of nutritional elements in the surrounding area of cells and therapeutic regimen [10].However, the mostly used models for the glioma case are the diffusive models that simulate the spatiotemporal change of glioma cell density, using partial differential equations.The combination of biomechanics and diffusive models is one of the latest advances in glioma modeling, for combining diffusion glioma cells and displacement of tissue caused by glioma growth [11].

THE DIFFUSION-REACTION MODEL
Unlike solid tumors, for which simple exponential or geometric expansion represents expansion of tumor volume, the glioma growth rate cannot be determined as the classical doubling rate [12], because gliomas can migrate and proliferate.In order to simulate glioma expansion, scientists have proposed the application of the diffusion-reaction equation, which is currently mostly used.The first that proposed a diffusive model for glioma growth was Murray in 1989 [13].() denotes the source factor, representing the glioma cell reproduction.
() denotes the treatment factor, representing the glioma cell loss due to treatment.This is zero, when no treatment is applied.
The initial state of the model, (, 0), is defined as the initial distribution of cancerous cells.
() ≡ () − () is the net proliferation rate.This equation was the basis of the later most important works for glioma modeling.The general procedure for glioma modelling, as derived by Murray's diffusive model, is presented in Figure 2.
The solution of the diffusion-reaction equation requires the application of numerical schemes, since there is no direct formula of its solution.The solution has to be approximated iteratively till time point of interest is reached.

Net Proliferation Rates for untreated gliomas
In 1995, Tracqui studied the evolution of cell concen- tration, by using two characteristic of tumour growth: proliferation and invasion [3].Tracqui proposed that the cells proliferate at exponential rate ρ, i.e.: So, Tracqui changed Equation 1, to where ρ denotes the proliferation rate of cells.Other proposed models [14], instead of geometrical rate, used either Verhulst law or Gompertz law where   is the maximum value that concentration can reach.Eq.2 has been mainly used for simulating untreated gliomas.

First Estimation of 𝑫𝑫 and 𝝆𝝆
From the very beginning of diffusive models, one of the key issues was the estimation of parameters that are being used.The diffusion coefficient  and the proliferation rate  are the basic parameters for the diffusive models.The first estimation of them is places in 1995 [1], when Silbergeld studied biological data and introduced two groups of glioma cells: the common ones and the resistant-to-first-chemotherapy ones.Parameter D was firstly defined either at  =  −   /, with the percentage of cells resistant to chemotherapy being at 8%, or at  =  −   / without resistant cells, while ρ was defined at  =  − /.

Resection Modeling
Getting this further, one of the next steps was to the model cancer evolution after ectomy [3,15].This was firstly simulated in 1993 by setting the concentration of the ectomized area equal to zero and, then, allowing the surrounding malignant cells proliferate and diffuse until the sphere reaches 6cm diameter.An example of ectomy, reproduced from [3], is given in Figure 3.

Low Grade Gliomas
Up to 1996, diffusive models studied high-grade gliomas due to their remarkably fast invasion.However, studying low-grade gliomas was important as well.Hence, in 1996, Woodward [16] suggested that speed of growth in low-grade tumours should be 10% of the respective one in high-grade gliomas, yielding satisfactory results.After some years, in 2003, Mandonnet [17] proposed that low-grade gliomas grew slowly, but linearly.This is mathematically derivable by Eq.3, because the expanding velocity of a population, which follows only the diffusion and growth laws of (3)

Brain Heterogeneity
Up to 2000, researchers didn't take brain anatomy into account.However, taking into account brain matter is of foremost importance, since it was observed that migration in white brain tissue is faster than in grey tissue [4].Thus, modeling of gliomas should take brain heterogeneity into account.Indeed, in 2000, Swanson innovated by introducing the problem and incorporating white and gray matter differentiation in the diffusion coefficient  of Eq.3.More specifically, the equation continued to hold, but for variable , i.e.: Diffusion coefficient  varies according to position, with D(x) =   or   , i.e. being constant for  in grey and white brain tissue respectively.Moreover, in order Swanson to apply the model, an brain atlas was required providing all the information about white and grey matter areas.Indeed, BrainWeb database [18] was available for extracting this information.An example of Swanson's simulation, reproduced from [4], is given in Figure 4.

Anisotropic Cell Migration
In 2005, Jbabdi et al. [19] introduced brain tissue anisotropy in diffusive modeling.As observed, glioma cell migration is facilitated along the directions of white matter fibers [20][21].This observation can be supported by the diffusion tensor Magnetic Resonance imaging (DT MRI) that gives a very good 3D reconstruction of white matter fibers.Jbabdi defined Eq.7, as where  is the diffusion tensor that describes cell diffusion rate at point , i.e. a 3 × 3 symmetric positive definite matrix that reflects local anisotropy.

Biomechanical Deformation
In 2005, Clatz et al. from INRIA [11] developed a model that simulates Jbadbi's model, but taking also into account the biomechanical deformations that occurs in brain due to tumour expansion.This model uses a predictor of the mass effect induced by both the tumor proliferation and infiltration.Figure 5 shows an example of a simulation with Clatz's model, with the deformation of tissues presented.The image is reproduced from [11].

Parameter estimation in modern models
According to [22], Figure 6 can be used as a guiding index for defining parameters  and  according to glioma grade, velocities of growth and /.Indeed, this log to log graph includes all parameter values found up to Jan. 2007 for both low-and high-grade gliomas.
Low-grade gliomas are sited in bottom left rectangle (LGG), for 2 mm/yr average velocities.Respectively, high grade gliomas (HCG) are positioned in the large rectangle, defined by / of 2 to 20 cm 2 and average velocities from 10mm/yr to 200mm/yr.On the left part gliomas with detectable mass, which can be cured with surgery, are placed.Continuing, for the case of heterogenous brain matter for high-grade gliomas, it is suggested in [4] that a typical value for  is =0.0012/day.This value for lowgrade gliomas can be defined at  = 0.00012 /day.Then by assuming that   = 5  ,   = 10 −2 mm 2 /day and   = 2 • 10 −3 mm 2 /day can be used [19].
It is noteworthy to present the following table, reproduced from [23], where there is a list of the proposed values of each parameter that the latest glioma models use.

CHEMOTHERAPY MODELING
Chemotherapy modeling is a quite blurred part of glioma models, which is mostly studied currently by researchers.Parameters that have to be taken into consideration Time interval (step)  1 day (CR) 6 hours (HFR) [29,30] Copyright © 2010 SciRes.JBiSE should be extracted by histopathological and biostatistical data of specific patients [31].Swanson has introduced a generalised net proliferation rate of chemotherapy as: where () is the temporal profile of the chemotherapy treatments, assuming a loss proportional to the strength or amount of therapy at a given time.Swanson sets () =  when the chemotherapy is being administered and () = 0 otherwise. is actually a measure of effectiveness of chemotherapy.Thus, in order therapy to be effective and size of tumor to decrease,  should be larger than , so that the net proliferation rate () is negative.This means the number of the dying cells is larger than the new born cells.
An example of chemotherapy application on real clinical MRI data, using ( 8) can be reproduced by [32].In this example, k was intentionally set to the high value of () = 0.024, so as the cancer to have shrinking effect.The data (18 MRI slices) has been acquired by Universität des Saarlandes Klinikum (Germany) within the scope of the Contracancrum project [33].The modeling results are given in Figure 7.The 3-dimensional representation of the initial and after-chemotherapy states are presented, accompanied by a sampled series of 4 MRI slices.

Mathematical Frameworks
The increasing interest of researchers on diffusive mod-els was remitted by the total absence on specific guidelines on how to design them mathematically.In 2009, Roniotis et al. [34] published a mathematical framework on designing diffusive models by using Finite Differences, for the needs of the ContraCancrum Project [33].The framework includes heterogeneous tissue, anisotropic migration of cells along white fibers and is 3D.Moreover, chemotherapy formalism of Eq.8 can also be incorporated in the model.
Different schemes of finite differences have been developed, namely forward Euler, backward Euler, Crank Nikolson and θ-methods.Moreover, the same model has been also developed with finite elements.The accuracy of the different schemes has been tested on simplified cases which suggested that the backward Euler scheme (Finite Differences) yields the best results.An example of applying this framework on real clinical data is displayed in Figure 8, reproduced from [34].

CHALLENGES FOR FUTURE MODELS
Simulations of applying therapy can be performed by adjusting the proliferation term according to glioma cell proliferation rates in the diffusion equation.Moreover, consecutive chemotherapy sessions have been simulated by works.However, the main issue governing these applications is the estimation of the most efficient model parameters.This requires study of pharmacokinetics and study of real chemotherapy sessions on patients.An ideal model will have the ability to suggest paths of invasion, unseen by the doctor and help the clinician predict how the tumor is going to behave after different cures.Personalized parameter estimation for treatment term, proliferation rate and initial state are very important for the model and are being studied by most researchers working on diffusive models globally.
Another important issue regarding modern diffusive models is the estimation of parameters of anisotropic diffusion.This is already incorporated in the models, but the methods used for extracting diffusion tensors have not been validated.This is mostly due to the use of DT-MRI, that is a modern technique and is hard to accompany all medical data.An estimation of tensors based on atlas could be studied.

CONCLUSIONS
Mathematical diffusive modeling for simulating glioma growth and invasion is a currently developing scientific field.An efficient 3-dimensional model of glioma growth will constitute a powerful tool for clinicians, since they could predict how glioma is going to develop in time.The most important part of this is the differentiation of the predicted outcome according to the different model input parameters.These parameters are changing with different therapy treatments.Thus the clinician can make an easier decision on which therapy scheme yields the best predictions.

Figure 1 .
Figure 1.Mathematical modeling for a general biological process.
Mandonnet et al. used clinical data reproduced from 27 patients to estimate that the average tumour velocity was 2 mm per year.

Table 1 .
[23]roduced from[23])-The parameters of the diffusive model as proposed in the last 2 decades.