cancer causes more than 13% of all deaths worldwide with an estimation of 12 million deaths in 2030 [ 1 , 2 ] . in united states , 2.5% of cancer deaths are caused by brain tumors , with grade iv astrocytoma , namely , glioblastoma multiforme ( gbm ) , accounting for 23.1% of them [ 46 ] . despite extensive research on gbm during the last decades , mortality has not changed significantly over the last years , with average life expectancy ranging 12 months after diagnosis and only 4% of treated patients being alive after 5 years [ 712 ] . unfortunately , the detection rates of the exact boundaries of gbm with common imaging techniques , such as magnetic resonance imaging ( mri ) , x - ray computed tomography ( ct ) , and positron emission tomography ( pet ) , are still poor [ 1214 ] . unlike solid tumors , for which simple exponential or geometric expansion could represent expansion of tumor volume , the gbm growth rate can not be determined as the classical doubling rate , because gbm does not form a solid and compact mass with cells invading the surrounding lesion . clearly , new mathematical formulations are necessary when studying this specific glioma case . since early 90s , there has been a vast amount of research towards simulating and formulating the mechanisms of gbm development , both in macroscopic and microscopic levels . microscopic models [ 1620 ] study the intracellular biological interactions in cell level , while macroscopic models [ 2127 ] study the tumor behavior , velocity , and mass deformation with using anatomical information derived from medical images . gompertzian and volumetric macroscopic models [ 2830 ] did not produce realistic clinical cancer growth representations as they lack taking invasion of tumor cells into account . the class of macroscopic models that has been widely used for simulating this diffusive behavior is diffusive modeling . these models make use of equations based on diffusion reaction schemes of glioma cell concentration , in order to couple the migration of gbm cells ( diffusion term ) and the proliferation of cells ( reaction / source term ) . the first diffusive gbm model was introduced in 1995 by tracqui   who used the diffusion reaction equation ( dre ) for simulating the spatiotemporal change of glioma cell concentration . the dre writes as follows : 
 ( 1)c(x , t)t = div(dc(x , t))+f(c(x , t ) ) , 
 where c(x , t ) is the tumor concentration in voxel x = ( x , y , z ) at time t , d is the diffusion coefficient ,  and div are the gradient and divergence operators respectively , and f(c ) is the net cell proliferation rate . later , swanson et al .   altered this equation by taking into account the high velocity of gbm in myelin sheaths by setting different diffusion coefficients in white and gray matter ( wm , gm ) . the new equation was changed to 
 ( 2)c(x , t)t = div(d(x)c(x , t))+f(c(x , t ) ) , 
 where d(x ) is the local diffusion coefficient , being either dg or dw for x being in gm and wm respectively . in that model , dw is assumed to be five times larger than dg for high grade glioma , according to the clinical observations of giese et al . [ 32 , 33 ] . the next change of the dre equation was proposed by price et al . who used t2- and diffusion tensor imaging ( dti ) mris to identify abnormalities caused by gbm [ 6 , 34 ] . therefore , jbabdi et al .   introduced diffusion tensors instead of gradient diffusion coefficients . the new advanced equations for simulating anisotropic growth of gbm in wm writes as follows : 
 ( 3)c(x , t)t = div(d(x)c(x , t))+f(c(x , t ) ) , 
 where d(x ) is the local diffusion tensor which describes the directional tumor cell diffusion , that is , a 3-by-3 symmetric positive definite matrix . continuing , clatz et al .   used biomechanics to simulate the deformation of brain structures caused by the development of gbm . in 2007 , hogea et al . introduced an advection term in ( 2 ) , with a two - way coupling with the underlying tissue elastic deformation . adapted the model to successive sessions of mris by using eikonal equations for coupling the reaction - diffusion dynamics with tumor evolution . the characteristics of gbm lead to aggressive treatment strategies that most often include surgery , irradiation and chemotherapy . the mathematical approaches that were used for simulating gbm growth inevitably raised the need for incorporating techniques for therapy simulation . including therapy parameters in gbm models could help the clinicians optimize therapy schemes , as they could predict the response of patients to different therapeutic plans . indeed , simulating therapy can be performed by adjusting the proliferation term f(c ) of the dre , which reflects the rate of proliferating cells ( b(c ) ) minus the dead cells due to treatment ( t(c ) ) . in order to simulate radiotherapy and chemotherapy swanson et al .   predicted survival time after resection by estimating the net rates of proliferation and diffusion , and their ratio from pre - treatment gadolinium - enhanced t1-weighted and t2-weighted mr images . simulated radiotherapy by introducing the radiobiology parameters  and  into t(t ) , which are interpreted biologically as repairable single and lethal double - strand breaks to the cell 's dna , respectively . lastly , some important works have been proposed for non diffusive models of glioma , taking into account the biological mechanisms involved in tumor and normal tissue [ 38 , 39 ] . the most commonly used formalisms of b(c ) are either the geometric law 
 ( 4a)b(c)=c 
 the logistic law : 
 ( 4b)b(c)=ccmccm , 
 or the gompertzian law : 
 ( 4c)b(c)=clncmc , where  is the geometrical growth parameter and cm is the maximum tumor cell concentration parameter . [ 40 , 41 ] used the proportions of white and gm for calculating the diffusion coefficient of the dre , rather than using discrete diffusion coefficients . the proportions of white and gm , as well as the diffusion tensor were extracted by the open sri-24 brain atlas , thus no dti processing was needed . this study extends our prior successes with in silico prediction of tumor growth to incorporate the effects of radiotherapy in real clinical datasets . we use the proportional model of wahl et al . and apply the linear quadratic model   for modeling radiotherapy . compared to previous publications dealing with tumor growth modeling , our approach includes several improvements : the use of proportional local diffusion coefficients extracted from brain atlases , instead of discrete diffusion coefficients , the application of isotropic diffusion in gm and anisotropic diffusion in wm on the proportional diffusive model , the application of the linear quadratic model on our previous model   for simulating radiotherapy , the simulation of radiotherapy in multiple fractions by using the linear quadratic model , the evaluation of radiotherapy simulation by using jaccard , dice , and volume similarity metrics . the use of proportional local diffusion coefficients extracted from brain atlases , instead of discrete diffusion coefficients , the application of isotropic diffusion in gm and anisotropic diffusion in wm on the proportional diffusive model , the application of the linear quadratic model on our previous model   for simulating radiotherapy , the simulation of radiotherapy in multiple fractions by using the linear quadratic model , the evaluation of radiotherapy simulation by using jaccard , dice , and volume similarity metrics . this study makes use of the proportional diffusion model presented in , which is an extended version of ( 3 ) with 
 ( 5)d(x)=d(x)w(x ) , 
 where w(x ) is a symmetric matrix that describes the anisotropy of cell diffusion along the 3-dimensional directions for each position x and d(x ) denotes the proportional diffusion coefficient . in this case , w(x ) writes as 
 ( 6)w(x)=[wx(x)000wy(x)000wz(x ) ] , 
 where wi(x )  , i = x , y , z is the directional diffusion weight , which denotes the contribution of the respective direction to the anisotropic migration of gbm cells in position x. thus , w denotes the contribution of each axis to the direction of white fibers , while d actually denotes the velocity of cells along this direction . d and  can be approximated by using the method of harpold et al . . in order to make a total use of the proportions pg of wm and pw of gm that can be provided by brain atlases , such as the sti-24 atlas [ 45 , 46 ] , d(x ) does not take two discrete values dg and dw . instead , the model uses a continuous pattern of coefficients by computing d(x ) as follows : 
 ( 7)d(x)=pg(x)dg+pw(x)dw . this paper extends this idea by improving ( 5 ) , so as to turn off anisotropic migration in gm tissue . in order to apply anisotropy only to wm tissues , the identity matrix i3 is used for changing ( 5 ) to the following formalism : 
 ( 8)d(x)=pg(x)dgi3  +  pw(x)dww(x ) the sri24 atlas [ 45 , 46 ] is an mri - based atlas of normal adult human brain anatomy , generated by registering images of 24 normal brains . this atlas provides the proportions pg of wm and pw of gm in each position of the brain , which can be registered to the clinical datasets of the patients . the data is provided in 155 slices with 240  240 pixels , for both stripped and unstripped skull . apart from the matter proportions , the atlas provides the dominant eigenvectors of the diffusion tensor . these eigenvectors have been extracted from dti , by computing the covariance matrices of the distribution of the 3d gaussian probability which simulates the diffusion of water . thus , they can be directly used for simulating the anisotropic migration of glioma cells along white matter fibers , as they represent the directions towards which the water diffusion extends mostly . by mapping clinical mri data to a brain matter atlas , it is possible to approximate the required tissue information / composition in the tumor area , since this is not possible to extract from the real mri . hence , using local diffusion coefficients on real mri images becomes possible , even if the lesion  hides  the underlying tissue . this paper investigates an extension of the preceding proportional glioma model to include the effects of the radiotherapy using the linear quadratic model . according to the model , the probability of cells surviving s following a single dose of radiation r(x , t ) was observed to follow this relationship : 
 ( 9)s(r)=exp(arr2 ) , 
 where linear  ( gy ) and quadratic  ( gy ) are the radiobiology parameters , which are interpreted biologically as repairable single and lethal double - strand breaks to the cell 's dna , respectively . in our case , where the clinical datasets have been extracted from patients with radiation in a number of fractions with the same dose , ( 9 ) turns to the following equation : 
 ( 10)s(r)=exp(arrr ) , 
 where r is the dose per fraction and r = nr is the total dose for n fractions . in general , fast - proliferating cells , like gbm , have a high tissue response /. in our experiment , we use a constant ratio / = 10 that is widely used in clinical applications for highly developing cancer and has been extracted from empirical dose response data . thus , if ( 4b ) is used for approximating cell proliferation , the overall proliferation term f(c ) for ( 3 ) turns out to be the following : 
 ( 11)f(c )  = { ( (1  ear(t)r(t)r(t)))ccmccm , ttherapy,ccmccm , ttherapy . the model uses t1- and t2-mris ( 255  255 pixel slices ) taken from six patients diagnosed with malignant glioblastoma multiforme and cured with radiotherapy . the datasets were provided by saarland university , germany , for the needs of the contracancrum project . for each dataset moreover , there are two or more sessions taken on different dates for tracking glioma development and evaluation simulation result . thus , they have to be firstly registered and interpolated to the atlas , before applying the dre equation . the optimized automatic registration 3d approach by mipav is used for registering the clinical datasets to the atlas . after registration , the datasets are delineated by expert radiologists from the university of saarland , by using doctoreye . figure 1(a ) presents the 99th slice of one patient dataset after registration and interpolation to the atlas . figure 1(b ) presents the 99th slice of the sri24 normal brain atlas . respectively , figures 1(d ) , 1(e ) , and 1(f ) present the 99th image of white matter proportion , gray matter proportion and the one dominant eigenvector . finally , figure 1(c ) presents the 99th slice of the mri set with the tumor area annotated by the clinician . in order to solve the dre , more specifically , the crank nikolson numerical scheme has been implemented for solving the dre in three dimensions . parameters dw ,  , and initial concentration of cells c0 , have been extracted from bibliography according to [ 41 , 44 ] . thus , we used  = 0.012 day , dw = 0.010  mm / day , dg = 0.002 for the solver of finite differences , we used a grid with dimensions x = y = z = 1  mm and a stable time step of t = 1  h. before applying the diffusive model , it is appropriate to estimate the values of the radiobiology parameters  and . the method used is that of rockne et al . , by projecting the temporal tumor curves for different values of a on the real datasets . estimating radius is ideal for spherical tumors however gbm has very complex shape thus , instead of tumor radii , the estimated total number of gbm cells is used for identifying the value of a that minimizes the error between simulated and actual target total gbm cells . figure 2 presents a set of simulated curves extracted for one specific patient after applying 7 different values of a. the real tumor total cells are marked on the curve for three different mri provided ( the three green points on day 1 , 15 , and 70 ) . the patient was treated with a total dose of 60  gy at 2  gy / fraction and 5 fractions / week for 6 weeks . each different curve is extracted by simulating radiotherapy for 0.003 increments of a , starting from zero ( no effect ) and reaching 0.021  gy . one can see that the optimal value for a , that minimizes the error from the total number of cells is a = 0.006  gy . the same process is followed for approximating the radiobiology parameters for the rest five patients . table 1 presents the patient specific parameters for all the patients . after defining all simulation parameters , figure 3 presents the simulated results for the 4th patient 21 days after the end of radiotherapy ( 100 days after diagnosis ) . the real dataset has been extracted on the same day of simulation . for this specific patient , radiotherapy caused shrinkage of the tumor , which can not be clearly seen on the simulation results of figure 3 . in order to have a better understanding of the shrinkage effect of radiotherapy , we use a line intersecting the point with the maximum concentration . then , concentration is estimated along this line and presented on a graph . this is applied both on the real initial mris ( day 1 ) and the simulated final result . figure 4(a ) shows the initial mri with the intersection line projected , while figure 4(b ) presents the graph of the concentration of gbm along these line . the red curve shows the initial estimated concentration along this line , while the blue curve shows the simulated concentration of radiotherapy along the same line . in order to have a better visualization of the shrinkage effect , we also provide figure 6 that depicts the change of the total number of cells in time ( for 100 days ) . we have made simulations on a series of six datasets in order to compare the simulated tumor with the final actual tumor , as annotated by the radiologists . for evaluation therefore , the annotated final tumor serves as a golden ground truth and the jaccard ( jc ) , dice ( ds ) and volume similarity ( vs ) metrics are adopted . the evaluation metrics are defined as : 
 ( 12)jc = tp(fp+tp+fn),ds=2tp(fp+2tp+fn),vs=1|fpfn|(fp+2tp+fn ) , 
 where tp ( true positive ) is the number of tumor voxels belonging to both ground truth and simulated result , fp ( false positive ) is the number of tumor voxels belonging to simulated result and not belonging to ground truth , and fn ( false negative ) is the number of tumor voxels belonging to ground truth but not belonging to simulated tumor . the average resulting jc , ds , and vs values for all patients are 95.21% , 98.09% and 99.62% . these values are 2.1% , 4.01% , and 1.89% higher than the respective results for discrete diffusion coefficients and 6.45% , 8.49% , and 4.07% than the uniform tumor growth model . despite the limited number of datasets , the evaluation results indicate that there might be a slight improvement in using the proposed model on glioblastoma multiforme . this is expected due to the several improvements introduced . by using the proportional diffusion coefficients , all available information is fully exploited for characterizing tissues , without needing to truncate tissue areas to either white or gray matter . moreover , the use of atlases contributes to defining the matter proportions on the areas underlying the tumor area . this could not be true without the use of atlas , as the tumor hides the tissue matter needed for defining the diffusion coefficients at the tumor area . also , by using anisotropy only on the portions of white matter , instead of all type of tissue , we better approximate the glioma migration along white fibers . radiotherapy is usually applied in fractions and the linear quadratic model has been adjusted to this , contrary to previous diffusive models which use the quadratic model for a single dose . by using the formalism of ( 10 )
applying diffusive models for simulating the spatiotemporal change of concentration of tumour cells is a modern application of predictive oncology . 
 diffusive models are used for modelling glioblastoma , the most aggressive type of glioma . 
 this paper presents the results of applying a linear quadratic model for simulating the effects of radiotherapy on an advanced diffusive glioma model . 
 this diffusive model takes into consideration the heterogeneous velocity of glioma in gray and white matter and the anisotropic migration of tumor cells , which is facilitated along white fibers . 
 this work uses normal brain atlases for extracting the proportions of white and gray matter and the diffusion tensors used for anisotropy . 
 the paper also presents the results of applying this glioma model on real clinical datasets .