A mathematical model of tumor growth and its response to single irradiation
 Yoichi Watanabe^{1}Email author,
 Erik L. Dahlman^{1},
 Kevin Z. Leder^{2} and
 Susanta K. Hui^{1}
DOI: 10.1186/s1297601600327
© Watanabe et al. 2016
Received: 1 December 2015
Accepted: 19 February 2016
Published: 27 February 2016
Abstract
Background
Mathematical modeling of biological processes is widely used to enhance quantitative understanding of biomedical phenomena. This quantitative knowledge can be applied in both clinical and experimental settings. Recently, many investigators began studying mathematical models of tumor response to radiation therapy. We developed a simple mathematical model to simulate the growth of tumor volume and its response to a single fraction of high dose irradiation. The modelling study may provide clinicians important insights on radiation therapy strategies through identification of biological factors significantly influencing the treatment effectiveness.
Methods
We made several key assumptions of the model. Tumor volume is composed of proliferating (or dividing) cancer cells and nondividing (or dead) cells. Tumor growth rate (or tumor volume doubling time) is proportional to the ratio of the volumes of tumor vasculature and the tumor. The vascular volume grows slower than the tumor by introducing the vascular growth retardation factor, θ. Upon irradiation, the proliferating cells gradually die over a fixed time period after irradiation. Dead cells are cleared away with cell clearance time. The model was applied to simulate pretreatment growth and posttreatment radiation response of rat rhabdomyosarcoma tumors and metastatic brain tumors of five patients who were treated with Gamma Knife stereotactic radiosurgery (GKSRS).
Results
By selecting appropriate model parameters, we showed the temporal variation of the tumors for both the rat experiment and the clinical GKSRS cases could be easily replicated by the simple model. Additionally, the application of our model to the GKSRS cases showed that the αvalue, which is an indicator of radiation sensitivity in the LQ model, and the value of θ could be predictors of the posttreatment volume change.
Conclusions
The proposed model was successful in representing both the animal experimental data and the clinically observed tumor volume changes. We showed that the model can be used to find the potential biological parameters, which may be able to predict the treatment outcome. However, there is a large statistical uncertainty of the result due to the small sample size. Therefore, a future clinical study with a larger number of patients is needed to confirm the finding.
Keywords
Mathematical model Tumor volume Gamma knife Radiation effectBackground
Mathematical modeling of biological processes is widely used to enhance quantitative understanding of biomedical phenomena. This quantitative knowledge can be applied in both clinical and experimental settings. One important application of modeling exercises is in the area of cancer biology [1, 2]. Many mathematical models have been developed to represent some aspects of cancer [3–6]. Those models vary from a simple model trying to simulate the growth of tumor volume to sophisticated models including many biologically important molecular processes [7, 8]. Macroscopic models can emulate clinically relevant phenomena for anticancer therapy. Because of the complexity of the biology, there is a need to link microscopic/molecular processes to those macroscopic ones that we observe in the clinic. Accordingly, there would be a need for a multiscale model rather than just a macroscopic one.
Recently, many investigators began studying mathematical models of tumor response to radiation therapy. The models used for those studies can be categorized into either stochastic or continuum model. Some did extensive studies of the growth and radiation response of the tumor using stochastic models including many biologically important processes at a cellular level, such as apoptosis, angiogenesis, and hypoxia [9–11]. Other groups used a continuum model. The tumor growth was represented as diffusion in a threedimensional space and the model was successfully applied to study the tumor growth and its response to radiation, mainly for brain tumors such as glioblastoma multiforme (GBM) and lowgrade glioma [12–14]. There were also several studies in which the nonstochastic modelling approach was chosen to investigate the tumor kinetics after irradiation [15–18]. In these studies, the tumor volume was represented as a volume in onedimensional space or the number of cancer cells, and the tumor volume was considered to consist of two components, i.e., proliferating cells and dead cells, after irradiation. Such a model was used to simulate the response of the tumor volume to radiotherapy.
In this study, we proposed a new simple onedimensional model for the time evolution of the tumor volume before and after a radiosurgical procedure. We assumed that the tumor volume is composed of actively dividing (or proliferating) cells and nondividing (or dead) cells after the irradiation. Such a twocomponent model was originally developed in the 1970s [19]. This model was previously used to study the radiation response of tumor during radiation therapy [15, 16, 20]. Our model is different from these “standard” models in at least two aspects. First, the slowdown process of tumor growth with increasing volume was modeled by assuming that the growth rate is a function of the ratio of vascular and tumor volumes and the growth rate of the vasculature is slower than that of the tumor itself. Secondly, our radiation response model was created by assuming that the cells do not die instantly, but survive for a few cell cycles. This is supported by some experimental results [21–24]. It is noteworthy that the possibility of one or more mitotic potential of damaged cells was recently considered in a mathematical model of lowgrade glioma treated with fractionated radiotherapy [13] although the authors eventually concluded that this effect is not significant for this treatment. In the current model, the radiationinduced cell kill rate was represented by introducing the probability that the cell stops dividing at the end of the cell cycle. The probability was related to the survival fraction of the linearquadratic (LQ) model since the probability is not readily available in the current literature.
To test the newly developed model, we exploited the model to reproduce the temporal volume change of rat rhabdomyosarcoma tumors for varying radiation dosages. In the 1960s to early 1970s, several British radiobiologists undertook extensive animal studies about the effect of ionizing radiation on tumor volume [25–28]. These studies provide useful data to test our model. Among those, the article of Barendsen et al. [25] contains experimental data showing the tumor volume before and after a single fraction dose of radiation. Additionally, the model was used to simulate the metastatic tumors of five patients who were treated for their metastatic brain cancers with Gamma Knife stereotactic radiosurgery (GKSRS) technique.
Our model is simple, but it includes two fundamental processes, the effect of the blood supply to the growth rate and the prolonged mitotic capacity of the cells damaged by radiation, which have not explicitly been included before together within a simple formalism. The applications of the model to clinical cases of GKSRS suggested some of the model parameters as a potential predictor of the treatment outcome.
The outline of this paper is as follows. In Section two, we present the mathematical model with our reasoning behind the selection of the proposed formulae. Numerical solution methods of the set of nonlinear ordinary differential equations are briefly discussed in the same section. We used two solution methods. One method employed manual iterative adjustment of the model parameters to fit the model with experimental and clinical data. For another method, we applied a simulated annealing technique to automatically estimate the parameters to achieve the best fit between the model and the measurement data. In section three, we present the sources of tumor volume data, to which our model was fitted. We used the data from the old literature of animal experiments and inhouse data of five patients treated by GKSRS. The results of the application of the model to those data will be shown in Section four. After discussion of the models and the results in Section five, we summarize this study in Section six.
Mathematical model
The proposed tumorvolume kinetic model
This model can be represented by a set of three equations of three variables: the volume of proliferating tumor, V_{T}(t), the volume composed of nondividing cells, V_{ND}(t), and the tumor growth rate λ(t), with four constants, a radiobiological parameter α, the initial tumor growth rate λ(0), the vascular growth retardation factor θ, and the cell clearance rate η_{cl}.
The model parameters α, λ(0), θ, and η_{cl} were obtained to achieve the best fitting of the temporal change of the total tumor volume, i.e., the sum of V_{T} and V_{ND}, between the model and the experimental/clinical data.
Derivation of governing equations
The function f is the growth rate of the malignant dividing cell and depends on both the total volume of actively dividing cells, V _{ T }, and the tumor vascular volume, V _{ v }. The function g is the transition rate of the dividing cells to the nondividing cells under the influence of radiation. Both f and g are a function of dose D. The second term in Equation (11) indicates that the nondividing dead cells are removed away from the original location with the rate of η _{ cl,} which is related to the cell clearance time T _{ cl } by 0.693/η _{ cl }.
In the current model, the dividing cells are the cells which can divide or proliferate regardless of its clonogenic capacity. Hence, the dividing cells include the cells eventually die due to radiation. Meanwhile, the nondividing cells are the cells which have stopped dividing and eventually disappear from the tumor volume. Note that we do not consider the cells in a quiescent state such as the cells in the reversible G0/G1 cell cycle arrest state [31].
Growth rate limited by blood volume
Here, λ is the growth rate, which is modulated by the dosedependent cell proliferation probability p. The function λ depends on the volumes of both the tumor and the vascular structure. Note that the tumor volume doubling time T _{ d } is defined as 0.693/λ in this study.
The constant θ is called the retardation factor of the vascular structure (or the vascular growth retardation factor). We considered that θ ≤ 1 to realize the saturation phenomena of tumor growth with increasing tumor volume. It is noteworthy that the use of Equations (15) and (16) along with Equation (10) with only the growth term in the righthand side of the equation leads to the wellknown Gompertzian solution of the tumor growth. See Additional file 1 for further discussion.
Radiationinduced cell death
Furthermore, V _{ ND }(t) = 0 for t ≤ t _{−} for single fraction treatment.
The number of surviving cells after irradiation of the current model can be the same as the standard model when τ _{ rad } is set approximately equal to 3T _{ m }. See Additional file 1 for details.
We assumed that the tumor growth rate is governed by Equation (15). During the period of τ _{ rad }, we set the right side of this equation to zero; thus assuming that the volume ratio of the vascular structure and the tumor is constant while the cells die due to the radiation effect.
Solution methods
A set of nonlinear ODEs, Equations (1, 2), (3, 4, 5), and (6, 7, 8), were numerically solved. The initial conditions were the volume at the first radiology exam for the tumor volume and zero for the volume of nondividing cells. The initial value of the tumor growth rate was obtained from the tumor doubling time, T_{d}(0). By introducing a characteristic time that was equal to the cell cycle time and the nominal volume of the dividing tumor cells, those equations were nondimensioned for the solution. The solution was obtained by the ODE solver of MATALB, ode15s, which was developed for solving stiff ODEs (The MathWorks Inc., Natick, MA). For the majority of cases in the current study, the parameters were iteratively adjusted so that the temporal change of the tumor volume obtained by the model fits the measured data the best through visual examination.
Here Γ is a set of feasible parameter values.
For the remainder of this study, we used that α/β = 10 days for tumor, the characteristic time T* = 1 day, and the observation time T_{m} = 10 days. The active radiation effecttime τ_{rad} is not an optimization parameter but is set to a value between 3 and 10 days. Currently, we do not know well how fast remaining alive cells grow while cells are dying because of irradiation. Hence, it was assumed that λ is constant during the time period of τ_{rad} and it keeps decreasing according to Equation (8) after that.
Quantitative evaluation of models
Here, k is the number of parameters used in a model. The model parameters which were optimized for the data fitting were α, θ, T_{d}, T_{cl}, and τ _{ rad }. It is noted, however, that only α and θ were varied and T_{d}, T_{cl}, and τ _{ rad } were constant for modelling the volume data of the rat experiments (see the next section), in which an experimently varied parameter was only the dose to the tumor.
Experimental and clinical data
Rat rhabdomyosarcoma experiments
The article of Barendsen et al. [25] contains experimental data showing the tumor volume before and after a single fraction dose of radiation. The dose levels considered were zero (or control), 1000, 2000, 3000, and 4000 rads (or cGy). For their study, rhabdomyosarcoma was grown in rats, which were treated either by 300 kV Xray or neutron beam. The data are reproduced as Figure 20.2A in the radiobiology textbook [29]. We used the 300 kV Xray data to test our model. It is noteworthy that the tumor volume was measured with a caliper by visually identifying the tumor lesion since no adequate imaging tool was readily available at that time.
Clinical examples: metastatic brain tumors treated by GKSRS
GKSRS is a common modality used to treat intracranial metastatic cancers. For GKSRS, we deliver a large single fraction of the dose to the target by 1.25 MeV gamma rays generated by radioactive cobalt60 (Co60) sources. The typical irradiation time is about a half hour though it depends on the size of the target and the age of the Co60 source. The typical dosage for malignant metastatic tumors increases as the tumor size decreases [37]. GKSRS patients have at least one pretreatment magnetic resonance imaging (MRI) scan. The treatment plan is created with MRI data taken on the day of the treatment. After the treatment, the patients have several followup MRI scans with the first one, usually, 2 months post GKSRS. Consequently, during the course of the treatment from the diagnosis to the termination of the followup, every GKSRS patient receives several MRI scans that use the Gadolinium (Gd) contrast enhanced T1weighted magnetic resonance imaging technique (CEMRI) with a spinecho pulse sequence. Therefore, we could collect clinical data showing the change of tumor volumes before and after GKSRS by examining the MRI data. We selected five cases among our patients. These cases were selected so that these can represent typical changes of the tumor volume clinically observed with GKSRS patients. With these clinical data, we investigated if our mathematical model could reproduce the characteristics of the time variation of the observed tumor volume by adjusting the model parameters.
There is a considerable uncertainty about what kind of substance is visualized with CEMRI. It is known that the vasculature of cranial tumors is leaky and the blood–brain barrier (BBB) breaks down [38]. Hence, the extracellular matrix around the cancer cells can absorb the contrast agent, resulting in a contrastenhanced volume. Meanwhile, the area around the dead nondividing cells produced by radiation may or may not enhance on CEMRI. For the current study, we simply assumed that the CEMRI volume is the sum of the dividing and nondividing cells.
Treatmentrelated parameters for clinical Gamma Knife stereotactic surgery cases
Clinical case number  

Parameters  Unit  1  2  3  4  5 
Primary cancer^{a}  NSCL  RCC  RCC  Testicular  Melanoma  
Prescription dose  Gy  22  16  18  16  20 
Prescription isodose line  %  50  50  50  50  80 
Maximum dose  Gy  44.0  32.0  36.0  32.0  25.0 
Mean dose  Gy  38.3  21.2  26.0  22.0  24.3 
Tumor volume at GKSRS  cm^3  0.26  5.50  0.69  5.98  0.87 
Tumor volume at end of followup  cm^3  0.073  0.21  0.20  4.82  23.60 
Total monitor duration  days  137  680  276  106  80 
Day of GKSRS from initial MRI  day  34  117  13  78  29 
Local control^{b}  CR  PR  PR  SD  PD 
One of the potential applications of a simple model such as ours is to discover a predictor(s) of treatment outcome. Hence, the predictive capability of the proposed model was examined for clinical GKSRS cases in the following way. Using the volume data from the model, first we calculated the change in the tumor volume, R_{40}, as the ratio of the volume on the 40th day from GKSRS to that at the time of GKSRS. Then, we evaluated the statistical correlation between R_{40} and five variables, i.e., the tumor volume doubling time at GKSRS, T_{d}(GKSRS), the αvalue, the θvalue, the tumor volume at GKSRS, V_{tumor}(GKSRS), and the mean target dose, D_{mean}. The correlation was analyzed by taking the nonlinear regression of the secondorder polynomial function between R_{40} and each of the variables. The Pvalue was estimated for each pair by using the linear regression analysis tool “fitlm”, which is available in MATLAB.
Results
Rat rhabdomyosarcoma experiments
Model parameters for rat rhabdomyosarcoma experiment data: the proposed model^{a}
Parameter\Case  Unit  Control  1000  2000  3000  4000 

α  Gy^{−1}  0.30  0.20  0.16  0.145  
Dose  Gy  10  20  30  40  
T_{d}(0)  days  1.35  1.35  1.35  1.35  1.35 
θ  0.72  0.72  0.74  0.795  0.838  
T_{cl}  days  4  4  4  4  
N  11  11  19  18  14  
MSD^{b}  0.01883  0.01099  0.00816  0.00832  0.03039  
AIC^{c}  −59.76  −66.00  −137.31  −128.23  −75.86  
S^{d}  1.00E + 00  9.83E01  2.84E01  5.38E02  4.51E03 
The model parameter optimization was also accomplished by using the simulated annealing method for the 3000 rad case. A solution was obtained by selecting the following feasible parameter domain Γ: α ∈(0.01, 0.2), θ∈(0, 1), T_{d} ∈(0.5, 15), and T_{cl} ∈(1, 50). We performed simulated annealing with 100000 steps to numerically solve the minimization problem given by Equation (31). The optimal solution was α = 0.1523 Gy^{−1}, θ = 0.7941, T_{d} = 1.5 days, and T_{cl} = 5.98 days. The residual error of this solution was 0.4327; whereas the residual error was 0.5219 when the parameters shown in Table 2 were used. Since the parameter values and the fitting quality are not significantly different and the numerical parameter optimization took much longer computing time, we decided to use the manual iterative method to obtain the optimized values of the model parameters in the rest of the current study.
Table 2 also shows the MSD values for five sets of the experimental data. We did the same modeling study with the standard model. The comparison between the two models is presented by comparing the MSD values in Additional file 2. The results clearly indicate the superiority of the new model over the standard model.
Clinical examples: metastatic brain tumors treated by GKSRS
Model parameters for radiation response: clinical GKSRS cases^{a}
Clinical case number  

Parameters  Unit  1  2  3  4  5 
α  1/Gy  0.09  0.19  0.19  0.1  0.05 
Tumor doubling time T _{ d } (0)  days  29.0  6.0  9.0  6.6  7.8 
Cell clearance time T _{ cl }  days  13.0  38.0  40.0  10.0  20.0 
Retardation factor θ  No dim.  0.62  0.77  0.53  0.78  0.99 
Mean dose  Gy  38.3  21.2  26.0  22.0  24.3 
Initial volume  cm^3  0.126  0.0065  0.271  0.031  0.101 
Irradiation day for GKSRS  day  34  117  22  78  29 
Tumor volume at GKSRS  cm^3  0.268  0.392  0.662  6.933  1.314 
Tumor doubling time at GKSRS  days  32.6  28.0  11.2  16.1  7.9 
Volume ratio, R_{40} ^{b}  0.16  0.57  0.58  0.69  10.0  
% Cell survival fraction  2.24  6.03  1.99  25.4  70.0 
 a)
The range of the α values is between 0.05 and 0.19. These values are on the smaller side of the commonly used range between 0.1 and 0.3 for tumors.
 b)
The tumor volume doubling time at the time of GKSRS, T_{d}(GKSRS), ranged from 7.9 days to 32.6 days.
 c)
The cell clearance time is between 10 and 40 days.
 d)
The vascular growth retardation factor θ is larger for the tumors that kept growing even after GKSRS. Particularly, the continued rapid growth after the treatment seen with case five required θ = 0.99, implying that the vasculature grew as rapidly as the tumor volume in this tumor.
Pvalue for predictors of response R_{40} ^{a}
Predictor  T_{d}(GKSRS)  α  θ  V_{tumor}(GKSRS)  D_{mean} ^{b} 

Pvalue: linear  0.269  0.246  0.0950  0.652  0.720 
Pvalue: 2^{nd} order  0.302  0.0264  0.0073  0.324  0.817 
Discussion
Cell growth model
It was demonstrated that the proposed model could reproduce the observed growth pattern regardless the initial tumor volumes for the clinical GKSRS cases. The initial volume doubling times of tumors, T_{d}(0), varied from 6.0 to 29.0 days. The volume doubling time of our model was not constant, but it increased with increasing tumor volume. At the time of GKSRS, the doubling time ranged from 7.9 to 32.6 days. Our own clinical data for the metastatic tumors treated by GKSRS show that the volume doubling time before GKSRS depends on the primary tumor type and the mean values are 42.5, 66.7, and 17.1 days for NSCL, renal cell, and melanoma, respectively [41]. The doubling time used for our model was within the range of those clinically observed data.
One of the major assumptions made for our model was that the growth rate of tumor volume is a function of the ratio of the volume of the tumor vascular structure to the tumor volume. For simplicity, we assumed that the rate is linearly proportional to the ratio of those volumes. We further assumed that the vascular structure grows with a growth rate proportional to the tumor growth rate. The proportionality constant is named as the vascular growth retardation factor, θ. This factor is considered constant and the θ value was between 0.53 and 0.99 in the clinical examples. The tumor growth rate, hence, decreases as the tumor size increases. It is noteworthy that the current hypotheses naturally lead to the wellknown Gompertzian growth of the tumor volume.
Radiation response model
The mechanism of cell death after irradiation is adequately summarized by an excellent review by B.G.Wouters in Chapter 3 of [23]. In response to radiation damage, cells lose the replicative capacity and die in genetically controlled mechanisms such as apoptosis (or programmed cell death), autophagy, senescence, and necrosis. These processes take place right at irradiation regardless of the cell cycle phase. However, it is considered that a major cell death mechanism in irradiated proliferating cells is a mitotic catastrophe, which leads to the cell death due to one of the cell death mechanisms mentioned above, but at the time of mitosis. Hence, this mechanism is also called reproductive or mitotic cell death. Which mechanism is a dominant death process depends on the type of cells, genetic status of the cells, and the microenvironment where the cells reside. Therefore, for the current study, we did not differentiate these mechanisms in the model.
We do not advocate the use of the standard radiation response model for two reasons. First, the survival fraction is a nondimensional value and the coefficient of the radiation term in the governing equation such as Equation (1) should have the unit of 1/Time. Second, the survival fraction was experimentally obtained by measuring the number of growing colonies in a petridish sometime after the irradiation. This time factor is not included in this expression. Hence, we proposed an alternative model, which used a function given by Equation (29), instead. The function g(D) is the cell killing probability that depends on the radiation dose and is inversely proportional to the counting time of the surviving colonies, T_{m}, which is the time equal to about ten cell cycles. In addition to the cell killing term, the dividing probability of cell is affected by radiation. This is represented by the probability p(D), or Equation (28). After the irradiation, the probabilities p(D) and q(D) remain constant for a certain period. This implies that the cell proliferation capacity is temporarily altered and the tumor cells continuously die during this period.
The probability of a cell to divide after irradiation, p, is represented as a function of dose as given by Equation (28). We adopted the LQ formula to relate p with the dose since we do not have alternative data for this relationship. However, for the parameter optimization, p itself can be used. In other words, the LQ model has little meaning in terms of modeling the actual radiation effects. This is one of the reasons that the α value of the clinical cases was smaller than the normal range.
Often, tumor volume increases after the treatment. The repair capability of the nondividing cells, i.e., sublethalcell repair mechanism [23, 29], helps to increase the tumor volume after irradiation. However, this repair mechanism is not sufficient to explain the volume growth because the repair takes place within a few hours after the irradiation [34]. The current study showed that the surviving tumor cells after the irradiation can explain the increase in the tumor volume. This phenomenon is often called as “accelerated repopulation”. There is another potential mechanism, which helps explain the tumor recurrence. It is known that tumor may be composed of proliferative cells and cells that are slowly dividing. The latter cells termed as “cancer stem cells (CSC)” exist for preserving selfrenewal and to protect from external insults such as chemotherapy and targeted therapies [42, 43]. It is worth pointing out here that there is strong evidence that CSC in glial tumors may grow after irradiation [44, 45]. Explicit modeling of this pathway was not attempted in the current study and the effect of CSC will be investigated in the future.
Our model was successfully fitted to the animal data. The difference among five data sets was only the magnitude of the dose that varied from zero to 40 Gy (or 4000 rad). These data contained many temporal data points before and after irradiation. The tumor volume was measured visually with the aid of biopsy; hence, the uncertainty caused by the MRIbased tumor delineation was reduced. Therefore, the application of our model to these data is an excellent method of model testing. The results demonstrated the soundness of the model more clearly than the results seen in the clinical cases.
Biological models and parameters
Our model has similarity to a twocomponent (or twocompartment) model that was developed in the 1970s [19]. The model assumes that the tumor volume is made of cells that are dividing (or proliferating) and cells that do not divide anymore because of fatal radiation damage. The nondividing cells that can be considered dead cells are transported away from the original location with some time constant, i.e., the clearance time. In the model, the tumor volume is a sum of those two types of cells. Hence, the rate of tumor volume shrinkage is governed by the clearance time. The model has been successfully applied to cervical and headandneck cancer cases [15–17]. We also adopted this cell clearance model in this study.
Recently, the twocomponent model was extended to a threecomponent model that includes a reversible state named as “cellcycle arrest” state and denoted by C [31] in addition to the proliferating cell (S) and dead cells in an irreversible state (P). Upon responding to external cellular stress, the cells in state P transit to state C. Some fraction of cells in C dies and moves to state S. However, the remaining fraction of cells goes back to P. The transition rates between the states depend on the stress level. The authors looked into the radiation as the stress to validate their hypothetical model. Although the experiments were done invitro using human fibroblast cells, their results suggest a possibility of another type of cell repair mechanism in addition to the sublethally damaged cell repair. Our model can be easily extended to test those hypothetical biological mechanisms in the future.
In this study, tumor volume was considered proportional to the number of cells by assuming the tumor cell size does not change during the entire course of treatment. We used the proportionality assumption to estimate the tumor volume change after irradiation by applying the cell survival data of classic radiobiology experiments, i.e., the LQ model. It should be pointed out that there is an experimental evidence which indicates a change in the cell size and the intercellular distance upon irradiation [28]; thus, leading to a breakdown of our assumption. Considering the strong dependence of such phenomena on the microenvironment and the tumor type, however, we used the current simplifying assumption. The effect of the cell size on the volume change will be studied in the future.
Estimated biological parameters
The rhabdomyosarcoma experiment was done by changing only the radiation dose. However, the best fitting of the model to the experimental data was achieved by varying the αvalue and the vascular growth retardation factor θ in addition to the dose. As the dose increased, the αvalue had to be decreased, whereas the θvalue had to increase as shown in Table 2. Note that the αvalue is related to the survival fraction of cells, S, after single 2 Gy exposure [34]. In fact, α = 0.145 and 0.3 with α/β = 10 lead to S = 71 and 49 %, respectively. In other words, a smaller α value implies less effective cell killing capability for the same dose. Thus, our modelling study suggests that the effectiveness of the irradiation decreased as the total dose increased. Note that this result contradicts with the standard LQ model, which would have a difficult time matching the experimental data with a fixed alpha since the dosages are so large. The increase in the θvalue was needed to reproduce the rapid growth of the tumor volume after irradiation. It is true that the θvalue before irradiation should not depend on the radiation dose. However, our model has only one θvalue. In the future model, it will be possible to assign different θvalue before and after irradiation for the model consistency.
To accurately model the response to the radiation for the clinical cases, we had to use α values (0.05 to 0.19 Gy^{−1}), which are on the smaller side of the commonly accepted range between 0.1 and 0.3 Gy^{−1}. The reason for the small α value is not clear, but it may be expected because the mean dose that was used for the modeling study may not be representative for GKSRS due to the highly nonuniform nature of the dose distributions. This warrants further investigation.
Delineation of tumor volume
In this study, we assumed that the volume visible as a contrastenhanced area on T1weighted MR images indicates the tumor. There is uncertainty in the volume delineated by this method, however. The Gadolinium contrast permeates into the exogenous and endogenous cellular space by the blood brain barrier disruption caused by the tumors [46]. It is a common assumption that the contrasted area contains both the proliferating cells and the dead cells. Hence, we used the sum of the dividing and nondividing cells as the tumor volume to be compared with the MRI data. It is not clear, however, if this disruption ceases to exist and BBB prevents the contrast material from transporting into the region when the cell dies.
Future directions
Our model contains eight biological parameters, i.e., α, α/β, T_{d} (tumor doubling time), T_{cl} (cell clearance time), θ (vascular growth retardation factor), τ_{rad} (effective time of radiation), T_{cc} (cell cycle time), and T_{m} (observation time or colony counting time). The radiation dose is known, but there is a possibility to use the dose as an unknown parameter because the physical dose is not uniform inside the tumor volume.
In this study, we showed that there was a significant correlation between the tumor volume change after GKSRS, R_{40}, and two model parameters: the αvalue and the θvalue. Since α indicates the sensitivity of the tumor to radiation, the result is easily predictable. But, the strong correlation between R_{40} and the θvalue is not obvious. This may imply that the tumor responds better to irradiation when the vascular structure grows much slower than the tumor volume. In other words, the shortage of the vasculature slows down or stops the regrowth of the tumor after irradiation. The current results were obtained from only five samples; hence, the confirmation of θ as a predictor of the GKSRS treatment outcome requires an extensive study with a much larger sample size. Such a future study will enable us to find a correlation between the local tumor control and the model parameters. If these parameters can be diagnostically obtained before the treatment, we should be able to design a patientspecific treatment by adjusting the treatment parameter, that is, the radiation dose to achieve a favorable treatment outcome.
There is experimental evidence showing that the radiation damage to the vascular structure may play a role in enhancing the cell killing capability of radiation for a single large dose of irradiation [47, 48]. For the current study, we assumed the vascular volume, V _{ v }, changes along with the tumor volume by assuming its growth rate proportional to the tumor growth rate. The proportional constant was given by the retardation factor, θ. Evidence of this phenomenon was observed with the GKSRS cases we studied as discussed in the above paragraph. To explicitly include the effect within our current modelling framework, we can easily add another equation for V _{ v } and include the radiation damage effect. However, this can be done only when more data on the magnitude of vascular damage as a function of the dose will become available.
Conclusions
We developed a simple mathematical model of tumor growth and its response to radiation by incorporating two key characteristics: (i) the tumor growth rate decreases as the tumor volume increases, and (ii) some radiationdamaged cells still keep dividing for a few more cell cycles after a single pulse of irradiation.
The current model was successful in mimicking both the animal experimental data and the clinically observed tumor volume changes. We showed that the volume changes of five tumors of Gamma Knife stereotactic radiosurgery patients could be fitted fairly well by using the proposed equations when appropriate model parameters were chosen. A correlation analysis indicated a strong relation between the postGKSRS tumor volume change and the α and θvalues.
Further refinement of the model which includes the radiationinduced vasculature damage is certainly desirable. Additionally, the model can be easily applied to a larger number of patients treated for GKSRS to find predictors/biomarkers of the treatment outcome. Such a study can confirm the importance of some of the modeling parameters introduced in this study, in particular, the θvalue, to predict better treatment outcome.
Methods
Computational calculations
All calculations were executed on a PC using MATLAB. The code is available from the first author.
Abbreviations
 BBB:

blood–brain barrier
 CEMRI:

contrast enhanced T1weighted magnetic resonance imaging
 CR:

complete response
 CSC:

cancer stem cells
 GKSRS:

Gamma Knife stereotactic radiosurgery
 LQ:

linearquadratic
 MRI:

magnetic resonance imaging
 MSD:

mean square of differences
 NSCL:

nonsmall cell lung
 ODE:

ordinary differential equation
 PD:

progressive disease
 PR:

partial response
 RCC:

renal cell carcinoma
 SD:

stable disease
Declarations
Acknowledgements
The part of this work was previously presented at the American Association of Physicists in Medicine (AAPM) annual meetings in the poster format, SUET08 and SUET309, in 2012 and 2013. The University of Minnesota Institutional Review Board approved the clinical data collection study (IRB #1103 M 97453).
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
Authors’ Affiliations
References
 Deisboeck TS, Zhang L, Yoon J, Costa J. In silico cancer modeling: is it ready for prime time? Nat Clin Pract Oncol. 2009;6(1):34–42. doi:10.1038/ncponc1237.PubMed CentralView ArticlePubMedGoogle Scholar
 Barillot E, Calzone L, Hupe P, Vert JP, Zinovyev A. Computational systems biology of cancer. Boca Raton: CRC Press; 2013.Google Scholar
 Kim Y, Magdalena AS, Othmer HG. A hybrid model for tumor spheroid growth in vitro I: theoreical development and early results. Math Models Methods Appl Sci. 2007;17:1773–98.View ArticleGoogle Scholar
 Cristini V, Li X, Lowengrub JS, Wise SM. Nonlinear simulations of solid tumor growth using a mixture model: invasion and branching. J Math Biol. 2009;58(4–5):723–63. doi:10.1007/s002850080215x.PubMed CentralView ArticlePubMedGoogle Scholar
 Deisboeck TS, Stamatakos GS. Editors. Multiscale cancer modeling. Chapman and hall/CRC mathematical and computational biology (book 34). Boca Rayton: CRC Press; 2010.Google Scholar
 Kim Y, Stolarska MA, Othmer HG. The role of the microenvironment in tumor growth and invasion. Prog Biophys Mol Biol. 2011;106(2):353–79. doi:10.1016/j.pbiomolbio.2011.06.006.PubMed CentralView ArticlePubMedGoogle Scholar
 Araujo RP, McElwain DL. A history of the study of solid tumour growth: the contribution of mathematical modelling. Bull Math Biol. 2004;66(5):1039–91. doi:10.1016/j.bulm.2003.11.002.View ArticlePubMedGoogle Scholar
 Powathil GG, Adamson DJA, Chaplain MAJ. Towards predicting the response of a solid tumour to chemotherapy and radiotherapy treatments: clinical insights from a computational model. PLoS Comput Biol. 2013;9(7):e1003120. doi:10.1371/journal.pcbi.1003120.PubMed CentralView ArticlePubMedGoogle Scholar
 Borkenstein K, Levegrün S, Peschke P. Modeling and computer simulations of tumor growth and tumor response to radiotherapy. Radiat Res. 2004;162(1):71–83.View ArticlePubMedGoogle Scholar
 Harting C, Peschke P, Borkenstein K, Karger CP. Singlecellbased computer simulation of the oxygendependent tumour response to irradiation. Phys Med Biol. 2007;52(16):4775–89. doi:10.1088/00319155/52/16/005.View ArticlePubMedGoogle Scholar
 Titz B, Jeraj R. An imagingbased tumour growth and treatment response model: investigating the effect of tumour oxygenation on radiation therapy response. Phys Med Biol. 2008;53(17):4471–88.PubMed CentralView ArticlePubMedGoogle Scholar
 Rockne R, Alvord E, Rockhill J, Swanson K. A mathematical model for brain tumor response to radiation therapy. J Math Biol. 2009;58(4):561–78. doi:10.1007/s0028500802196.PubMed CentralView ArticlePubMedGoogle Scholar
 PerezGarcia VM, Bogdanska M, MartinezGonzalez A, BelmonteBeitia J, Schucht P, PerezRomasanta LA. Delay effects in the response of lowgrade gliomas to radiotherapy: a mathematical model and its therapeutical implications. Math Med Biol. 2014. doi:10.1093/imammb/dqu009
 Nawrocki S, ZubikKowal B. Clinical study and numerical simulation of brain cancer dynamics under radiotherapy. Communications in Nonlinear Science and Numerical Simulation. 2014(0). doi:10.1016/j.cnsns.2014.08.001
 Lim K, Chan P, Dinniwell R, Fyles A, Haider M, Cho YB, et al. Cervical cancer regression measured using weekly magnetic resonance imaging during fractionated radiotherapy: radiobiologic modeling and correlation with tumor hypoxia. Int J Radiat Oncol Biol Phys. 2008;70(1):126–33. doi:10.1016/j.ijrobp.2007.06.033.View ArticlePubMedGoogle Scholar
 Huang Z, Mayr NA, Yuh WT, Lo SS, Montebello JF, Grecula JC, et al. Predicting outcomes in cervical cancer: a kinetic model of tumor regression during radiation therapy. Cancer Res. 2010;70(2):463–70. doi:10.1158/00085472.CAN092501.PubMed CentralView ArticlePubMedGoogle Scholar
 Chvetsov AV. Tumor response parameters for head and neck cancer derived from tumorvolume variation during radiation therapy. Med Phys. 2013;40(3):034101. doi:10.1118/1.4789632.View ArticlePubMedGoogle Scholar
 Zhong H, Chetty I. A note on modeling of tumor regression for estimation of radiobiological parameters. Med Phys. 2014;41(8):081702. doi:10.1118/1.4884019.PubMed CentralView ArticlePubMedGoogle Scholar
 Okumura Y, Ueda T, Mori T, Kitabatake T. Kinetic analysis of tumor regression during the course of radiotherapy. Struct Bond (Berlin). 1977;153(1):35–9.Google Scholar
 Chvetsov AV, Dong L, Palta JR, Amdur RJ. Tumorvolume simulation during radiotherapy for headandneck cancer using a fourlevel cell population model. Int J Radiat Oncol Biol Phys. 2009;75(2):595–602. doi:10.1016/j.ijrobp.2009.04.007.View ArticlePubMedGoogle Scholar
 Curtis SB, Barendsen GW, Hermens AF. Cell kinetic model of tumour growth and regression for a rhabdomyosarcoma in the rat: undisturbed growth and radiation response to large single doses. Eur J Cancer. 1973;9(2):81–7.View ArticlePubMedGoogle Scholar
 Forrester HB, Vidair CA, Albright N, Ling CC, Dewey WC. Using computerized video time lapse for quantifying cell death of Xirradiated rat embryo cells transfected with cmyc or cHaras. Cancer Res. 1999;59(4):931–9.PubMedGoogle Scholar
 Joiner M, van der Kogel A, editors. Basic clinical radiobiology. 4th ed. London: Hodder Arnold; 2009.Google Scholar
 Sakashita T, Hamada N, Kawaguchi I, Ouchi NB, Hara T, Kobayashi Y, et al. A framework for analysis of abortive colony size distributions using a model of branching processes in irradiated normal human fibroblasts. PLoS One. 2013;8(7):e70291. doi:10.1371/journal.pone.0070291.PubMed CentralView ArticlePubMedGoogle Scholar
 Barendsen GW, Broerse JJ. Experimental radiotherapy of a rat rhabdomyosarcoma with 15 MeV neutrons and 300 kV xrays. I. Effects of single exposures. Eur J Cancer. 1969;5(4):373–91.View ArticlePubMedGoogle Scholar
 Hermens AF, Barendsen GW. Changes of cell proliferation characteristics in a rat rhabdomyosarcoma before and after xirradiation. Eur J Cancer. 1969;5(2):173–89.View ArticlePubMedGoogle Scholar
 Thompson LH, Suit HD. Proliferation kinetics of xirradiated mouse L cells studied WITH TIMElapse photography. II. Int J Radiat Biol Relat Stud Phys Chem Med. 1969;15(4):347–62.View ArticlePubMedGoogle Scholar
 Tannock I, Howes A. The response of viable tumor cords to a single dose of radiation. Radiat Res. 1973;55(3):477–86.View ArticlePubMedGoogle Scholar
 Hall EJ, Giaccia AJ. Radiobiology for the radiologist. 7th ed. Philadelphia: Lippincott Williams and Wilkins; 2011.Google Scholar
 Ribba B, Kaloshi G, Peyre M, Ricard D, Calvez V, Tod M, et al. A tumor growth inhibition model for lowgrade glioma treated with chemotherapy or radiotherapy. Clin Cancer Res. 2012;18(18):5071–80. doi:10.1158/10780432.CCR120084.View ArticlePubMedGoogle Scholar
 Schäuble S, Klement K, Marthandan S, Münch S, Heiland I, Schuster S, et al. Quantitative model of cell cycle arrest and cellular senescence in primary human fibroblasts. PLoS One. 2012;7(8):e42150. doi:10.1371/journal.pone.0042150.PubMed CentralView ArticlePubMedGoogle Scholar
 Rockne R, Rockhill JK, Mrugala M, Spence AM, Kalet I, Hendrickson K, et al. Predicting the efficacy of radiotherapy in individual glioblastoma patients in vivo: a mathematical modeling approach. Phys Med Biol. 2010;55(12):3271–85.PubMed CentralView ArticlePubMedGoogle Scholar
 Puck TT, Marcus PI. Action of xrays on mammalian cells. J Exp Med. 1956;103(5):653–66. doi:10.1084/jem.103.5.653.PubMed CentralView ArticlePubMedGoogle Scholar
 Dale RG, Jones B, editors. Radiobiological modelling in radiation oncology. Oxfordshire: The British Institute of Radiology; 2007.Google Scholar
 Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical recipes: the Art of scientific computing. 3rd ed. New York: Cambridge University Pres; 2007.Google Scholar
 Brockwell PJ, Davis RA. Time series: theory and methods. 2nd ed. New York: Springer; 1991.View ArticleGoogle Scholar
 Shaw E, Scott C, Souhami L, Dinapoli R, Kline R, Loeffler J, et al. Single dose radiosurgical treatment of recurrent previously irradiated primary brain tumors and brain metastases: final report of RTOG protocol 90–05. Int J Radiat Oncol Biol Phys. 2000;47(2):291–8.View ArticlePubMedGoogle Scholar
 Greenwood J. Mechanisms of blood–brain barrier breakdown. Neuroradiology. 1991;33(2):95–100. doi:10.1007/BF00588242.View ArticlePubMedGoogle Scholar
 Eisenhauer EA, Therasse P, Bogaerts J, Schwartz LH, Sargent D, et al. New response evaluation criteria in solid tumours: revised RECIST guideline (version 1.1). Eur J Cancer. 2009;45(2):228–47.View ArticlePubMedGoogle Scholar
 Lin HY, Watanabe Y, Cho LC, Yuan J, Hunt MA, Sperduto PW, et al. Gamma knife stereotactic radiosurgery for renal cell carcinoma and melanoma brain metastases—comparison of dose response. J Radiosurg. 2013;2(3):193–207.Google Scholar
 Dalhman E, Watanabe Y. How fast do metastatic tumors grow in brain? 16th International Leksell Gamma Knife Society Meeting; 25–29 March. Sydney: Leksell Gamma Knife Society; 2012.Google Scholar
 Dick JE. Stem cell concepts renew cancer research. Blood. 2008;112(13):4793–807. doi:10.1182/blood200808077941.View ArticlePubMedGoogle Scholar
 Li L, Bhatia R. Stem cell quiescence. Clin Cancer Res. 2011;17(15):4936–41. doi:10.1158/10780432.CCR101499.PubMed CentralView ArticlePubMedGoogle Scholar
 Leder K, Holland EC, Michor F. The therapeutic implications of plasticity of the cancer stem cell phenotype. PLoS One. 2010;5(12):e14366. doi:10.1371/journal.pone.0014366.PubMed CentralView ArticlePubMedGoogle Scholar
 Leder K, Pitter K, Laplant Q, Hambardzumyan D, Ross BD, Chan TA, et al. Mathematical modeling of PDGFdriven glioblastoma reveals optimized radiation dosing schedules. Cell. 2014;156(3):603–16. doi:10.1016/j.cell.2013.12.029.PubMed CentralView ArticlePubMedGoogle Scholar
 Roberts TP, Chuang N, Roberts HC. Neuroimaging: do we really need new contrast agents for MRI? Eur J Radiol. 2000;34(3):166–78.View ArticlePubMedGoogle Scholar
 GarciaBarros M, Paris F, CordonCardo C, Lyden D, Rafii S, HaimovitzFriedman A, et al. Tumor response to radiotherapy regulated by endothelial cell apoptosis. Science. 2003;300(5622):1155–9. doi:10.1126/science.1082504.View ArticlePubMedGoogle Scholar
 Song CW, Park H, Griffin RJ, Levitt SH. Radiobiology of stererotactic radiosurgery and stereotarctic body radiation therapy. In: Levitt SH, editor. Technical basis of radiation therapy, medical radiology, radiation oncology. Berlin: Springer; 2012. p. 51–61.Google Scholar