Computer Implementation of a New Therapeutic Model for GBM Tumor

Modeling the tumor behavior in the host organ as function of time and radiation dose has been a major study in the previous decades. Here the effort in estimation of cancerous and normal cell proliferation and growth in glioblastoma multiform (GBM) tumor is presented. This paper introduces a new mathematical model in the form of differential equation of tumor growth. The model contains dose delivery amount in the treatment scheme as an input term. It also can be utilized to optimize the treatment process in order to increase the patient survival period. Gene expression programming (GEP) as a new concept is used for estimating this model. The LQ model has also been applied to GEP as an initial value, causing acceleration and improvement of the algorithm estimation. The model shows the number of the tumor and normal brain cells during the treatment process using the status of normal and cancerous cells in the initiation of treatment, the timing and amount of dose delivery to the patient, and a coefficient that describes the brain condition. A critical level is defined for normal cell when the patient's death occurs. In the end the model has been verified by clinical data obtained from previous accepted formulae and some of our experimental resources. The proposed model helps to predict tumor growth during treatment process in which further treatment processes can be controlled.


Introduction
Human brain is composed of cells called glios that proliferate to repair its injured part. This reproduction is totally under control of the brain. However, if any of these cells, without brain control, starts doubling for any reason, it will generate an unwanted tumor. GBM or glioblastoma multiform is the most malignant and aggressive type of high-grade gliomas which incur major clinical problems and mostly leads to death.
There are two types of brain cancer. Primary tumors are those that originate from the brain itself. The second type originates from metastasis in other parts of patient's body. The primary tumors may become cancerous or not; however, it needs room to grow; therefore, some types of injuries occur. Despite many advancements in surgery, radiotherapy, and other treatment methods, no effective treatment have been proposed in the previous decades and most of patients die of this type of cancer. Thus, primary brain tumors bring high mortality rate, which is attributed to short doubling time of the tumor. The treatment of this tumor encounters a series of problems. The tumor does not respond to schemes of treatment as expected. Another factor is related to the treatment plans that are not that much impressive and finally in major number of cases mathematical modeling is not used to provide an optimal dose delivery to the injured tissue.
Radiotherapy, being one of the main clinical solutions for cancer treatment, is technically the second remarkable method after surgery. Widespread application of this method is due to several factors like high precision, fast response, high-tech apparatus, and many decades of knowledge and experience. However, there are some disadvantages. Sometimes normal tissues are inevitably hit by radiation beams during dose delivery to the cancerous tissues. In order to 2 Computational and Mathematical Methods in Medicine gain a better view of injuries to normal and cancerous cells, a mathematical model is required. In the present paper, a new model of GBM tumor has been generated that does not suffer from deficiencies existing in other models.
GBM or glioblastoma multiform is the most malignant and aggressive type of high-grade gliomas which follow major clinical problems and mostly in 10 to 12 months leads to death [1]. Harpold et al. introduced [2] a model with radiotherapy input which was improved at the same year. Consequently a continuum mathematical model about the invasive cell treatment was developed [3]. A remarkable masterpiece of work was created for the first time in 2007 [4] which followed a two-state space equation without considering the input term of radiotherapy or chemotherapy. Tanaka et al. [5] worked out a hybrid model that showed that the proliferation of cancerous cells is dependent on the radius of the tumor. In 2010 [6] a model with the input term was written that used the radiobiological aspects. Based on this achievement, Barazzuol et al. [7] pushed the model further, considering the radiotherapy and chemotherapy but in a complicated way. Specifications of this model are investigated in the next section.

Model Development
Detailed information about properties and disadvantageous of abovementioned models is investigated.
Tanaka and his colleagues [5] suggested a hybrid compartment-continuum-discrete (CCD) model in order to simulate the proliferation of gliomas and the cell invasion. It is described by the equation below: One deficiency of this model is lack of state space equations and thus control over the system (tumor growth) is not easy. The other negative point is that the time parameter is missing and the future status of the cells cannot be estimated. Harpold et al. have written a review over previous models and tried to cast the proliferation of gliomas in a model. Their final result is model appearing in (2). Outcomes have slight difference with the aforementioned model [2]. Consider Harpold and his team continue their work and in 2007 their new model that has the net proliferation of gliomas cells [2] is where / is the rate of change of gliomas cell concentration and ( ∇ 2 ) is the net invasion of glioma cells and (1−( / )) is the net proliferation of gliomas cells. It is observed that the model contains time parameter and estimation of the amount of cancer cells is accessible but the disadvantage is the difference between mathematical results and experimental data. Besides, the status of normal cells is not simultaneously investigated and the death moment is not computable.
Stein et al. [3], in the same year in a study with the aim of defining the invasive cells treatment in glioblastoma, used a mathematical model called "continuum mathematical model" which is as follows: Its disadvantages include using too many parameters that have negligible effects on the final result plus the deficiencies of former models. Another work we want to mention is about what Rockne and his team did in 2008 [8]. They used the Swanson model and added the radiobiological effect by using a ( , ) term to the basic model of Swanson. They proposed the below model in a way that the LQ model of radiobiology is a part of their work: where / is the rate of change of gliomas cell concentration, ∇ ⋅ ( ( )∇ ) is the net dispersal of gliomas cells, is the net proliferation of gliomas cells, and ( , ) is the loss due to therapy.
Later in 2010 in a study with 9 patients [6] in whom the tumor had been diagnosed soon enough and underwent radiotherapy treatment, their model was improved on the basis of radiobiology of each patient for the increase of GBM cells. This model contains the LQ model of radiology in a more effective manner. This model responds to the treatment and for the first time the GBM model was combined with dose delivery input, which makes further investigation on the type of treatment schemes possible. Consider The first three terms are the same as above and the last term ( , ) is modified as follows: To obtain a general view of disadvantages of existing models, we can put forward points as using too many parameters that Computational and Mathematical Methods in Medicine 3 have negligible effects on its final result. The models' mathematics is not that much precise and is based on experimental data. The time parameter is missing and so estimating the future status of the cells is not available. The status of normal cells is not simultaneously investigated and the decease time of the patient is not computable. Most of papers here just have many statistics and they lack computation. In our model we provide new inputs in which the above deficiencies will be resolved. There are also novel issues in the proposed model.

Methodology
In the previous section, the deficiencies of former model have been declared. The present model consists of new parts that compensate most of the disadvantages of the other existing works.
What Rockne et al. [6] and Barazzuol et al. [7] did in 2010 is the basis of our idea. Their work is really effective and describes the status of normal and cancerous cell through time, but still its equations have just initial condition and lack the specific input. Here is their formula in two equations is the number of normal cells and is the cancerous ones. The other parameters are defined in the nomenclature part of this paper. Obviously, one can see that the states of and are not related to radiation dose input.
Although many of the other model deficiencies, like dependency on time, the stats of normal cells, and state space modeling, have been considered in this model, still there are two main drawbacks in this model. It has not been associated with an input term of radiotherapy. The effectiveness of dose for normal and cancerous cells is not defined. Barazzuol and other authors in 2010 [7] mentioned that according to the LQ model (see Appendix section) each cell (normal or cancerous) is diminished facing the treatment radiation and the main goal of this model is to open a new pathway for control issues. Since the Kirkby model lacks an input term, the amount of the interval of the treatment scheme for the GBM patients is not under control. If the input dose of radiation through this novel model is mathematically defined instead of equal nonprecisely calculated pulses of radiation a highly variable source of the decision making will be available for technical scientists to investigate ideas of their own about generating new plans in GBM treatment schemes.
There are four series of experimental data set that helped us in the GEP algorithm to generate the proposed model. In each set of data a specific portion of dose is delivered to the patients and various results have been accessed and are shown in Tables 2 and 3. It is easily seen that hyperfractionation (dose delivery is fractionated in more than once a day during treatment) affects the proliferation of normal and cancer cells in a way that death occurs after a longer period which prolongs lifetime of patient. In the first set of data as can be seen in Table 2 the whole treatment dose is delivered consistently. The second table shows a patient data that received the treatment in four parts of hyperfractionation, which obviously lead to a longer lifetime for the patient. Gene expression programming (GEP) similar to genetic algorithms (GA) and genetic programming (GP) is a gene based algorithm. Firstly by using populations of individuals, it selects them considering their fitness and introduces genetic variation using one or more genetic operators [9].
GEP uses chromosome's character linearly organized in a head and a tail made of genes. By means of mutation, the chromosome functions as a genome and is exposed to modification, transposition, root transposition, gene transposition, gene recombination, and one-and two-point recombination [10].
In the beginning (7) should be combined with LQ model (see Appendix section) in order to obtain a new differential equation which indicates the number of cancerous and normal cells of brain in the injured area. The equation should manifest the dose amount and the interval of each portion of dose delivery simultaneously. As it is said before (8) with a slight change is as follows: In the first step the above equation must be digitized and converted to a discrete form like this Based on (2) and the LQ equation (see Appendix section), the above information is given to the regression equation of GEP as the input or initial function. The objective of the present work is the discovery of a symbolic expression that satisfies a set of fitness cases. First, the set of functions and the set of terminals must be chosen [10].
Based on the above statements the input functions of the GEP algorithm are defined as follows. The four main operations as +, −, * , and / and ∧, also exponential functions like and ln( ), the triangular functions sin( ), cos( ), and tan( ) and hyperbolic functions sinh( ), cosh( ), and  In the algorithm utilized in our work, there are 18 operations (mathematical functions) and 5 operands (variables), summing up to 23. The nearest number to 23 in a binary format GEP is 32 which is 25. The length of the matrix based on 2 is 5 × 31 or 155 given the reason that 31 is the summation of the whole set. The initial number is 512, so the total dimension of the data is 155 × 512. A map for converting a digit more than 23 to a number between 1 to 23 is required here, since the operations such as mutation and crossover provide a digit more than the 23, the mutation probability is 0.1, and the crossover is supposed to be 13. The model is then capable of estimating the future status of the patient's treatment, if the proposed data, (8), and LQ model according to clinical dataset are combined by the GEP algorithm. The least square method (LSM) is then applied to choose the best fitting model gained from the GEP. The chosen sample is used for the regression and the obtained function is tested. There have been 100 sets of sample data divided into 4 groups. The first category has its dose delivery once a day.
The second one has its dose delivery twice a day and subsequently the third one three times a day and the last one is three times a day but with nonequal amount of dose.
Among each of 25 samples data, the first 20 have been used for training and the remaining 5 for testing. Since in each set of 25 samples data 7 to 10 of them treat totally different of the proposed model, data mining and PCA algorithm are applied to eliminate them from the training and test process. At the end the calculated model by using (7) and LQ model in their best form is suggested as follows: If the relations achieved by GEP are converted to a differential equation the following is generated: The formula parameters are as follows in Table 1.
In one portion of dose delivery each day the amount of dose was 1 Gy which was reduced to 0.5 Gy in the next categories.

Model Verification
In order to verify this model, a set of clinical data is required. Previous articles and our experimental tests provide this most important part of the present work. The model is fitted to attained data and the Kaplan-Meier [9] survival curves have been used. The weighting was adjusted heuristically to guide the optimization to fit correctly the tail of the curve representing the long-term survivors.

Results
Tables 2 and 3 are the tables of dose delivery for two separate patients as a sample of the clinical data. Table 2 is related to a patient that received 1 Gy of dose once a day. Table 3 is the data of another patient who had 0.5 Gy of radiation dose 4 times a day in different scheme of treatment.
All the results are depicted in Figures 1, 2, 3, 4, and 5. Each figure compares the error between the model response and the clinical data.

Conclusion
It is easily seen that this new model can efficiently empower the radiation oncologist and medical physicist control over the treatment procedure in comparison with the traditional models. Since new features according to the clinical and experimental data have been added to the proposed model, better and more reliable results have been generated. Here are some suggestions.  Increase in the number of clinical data can lead to more precise results. This investigation experienced around 95.2 percent of accuracy while fitting the theoretical data on the clinical ones. This value is obtained by using the average of error between each point of clinical and mathematical data.
We use a few variable parameters in this model, but one can obtain better formulae using a set of other parameters which are effective in cancer growth and modification.
Finally, the best further work is to generate a new intelligent algorithm to have new schemes of radiotherapy in   A proportional coefficient doseEff: The effect of dose on normal and cancer cells.