A mathematical model of tumor growth and its response to single irradiation
© Watanabe et al. 2016
Received: 1 December 2015
Accepted: 19 February 2016
Published: 27 February 2016
Mathematical modeling of biological processes is widely used to enhance quantitative understanding of bio-medical 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.
We made several key assumptions of the model. Tumor volume is composed of proliferating (or dividing) cancer cells and non-dividing (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 pre-treatment growth and post-treatment radiation response of rat rhabdomyosarcoma tumors and metastatic brain tumors of five patients who were treated with Gamma Knife stereotactic radiosurgery (GKSRS).
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 post-treatment volume change.
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.
KeywordsMathematical model Tumor volume Gamma knife Radiation effect
Mathematical modeling of biological processes is widely used to enhance quantitative understanding of bio-medical 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 anti-cancer 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 three-dimensional 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 low-grade glioma [12–14]. There were also several studies in which the non-stochastic 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 one-dimensional 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 one-dimensional 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 non-dividing (or dead) cells after the irradiation. Such a two-component model was originally developed in the 1970s . 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 slow-down 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 low-grade glioma treated with fractionated radiotherapy  although the authors eventually concluded that this effect is not significant for this treatment. In the current model, the radiation-induced 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 linear-quadratic (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.  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 non-linear 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 in-house 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.
The proposed tumor-volume kinetic model
This model can be represented by a set of three equations of three variables: the volume of proliferating tumor, VT(t), the volume composed of non-dividing cells, VND(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 VT and VND, 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 non-dividing 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 non-dividing 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 non-dividing 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 .
Growth rate limited by blood volume
Here, λ is the growth rate, which is modulated by the dose-dependent 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 right-hand side of the equation leads to the well-known Gompertzian solution of the tumor growth. See Additional file 1 for further discussion.
Radiation-induced 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.
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 non-dividing cells. The initial value of the tumor growth rate was obtained from the tumor doubling time, Td(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 non-dimensioned 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 Tm = 10 days. The active radiation effect-time τ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 α, θ, Td, Tcl, and τ rad . It is noted, however, that only α and θ were varied and Td, Tcl, 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.  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 X-ray or neutron beam. The data are reproduced as Figure 20.2A in the radiobiology textbook . We used the 300 kV X-ray 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 cobalt-60 (Co-60) sources. The typical irradiation time is about a half hour though it depends on the size of the target and the age of the Co-60 source. The typical dosage for malignant metastatic tumors increases as the tumor size decreases . GKSRS patients have at least one pre-treatment 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 follow-up 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 follow-up, every GKSRS patient receives several MRI scans that use the Gadolinium (Gd) -contrast enhanced T1-weighted magnetic resonance imaging technique (CEMRI) with a spin-echo 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 . Hence, the extracellular matrix around the cancer cells can absorb the contrast agent, resulting in a contrast-enhanced volume. Meanwhile, the area around the dead non-dividing 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 non-dividing cells.
Treatment-related parameters for clinical Gamma Knife stereotactic surgery cases
Clinical case number
Prescription isodose line
Tumor volume at GKSRS
Tumor volume at end of follow-up
Total monitor duration
Day of GKSRS from initial MRI
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, R40, 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 R40 and five variables, i.e., the tumor volume doubling time at GKSRS, Td(GKSRS), the α-value, the θ-value, the tumor volume at GKSRS, Vtumor(GKSRS), and the mean target dose, Dmean. The correlation was analyzed by taking the non-linear regression of the second-order polynomial function between R40 and each of the variables. The P-value was estimated for each pair by using the linear regression analysis tool “fitlm”, which is available in MATLAB.
Rat rhabdomyosarcoma experiments
Model parameters for rat rhabdomyosarcoma experiment data: the proposed modela
1.00E + 00
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), Td ∈(0.5, 15), and Tcl ∈(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, Td = 1.5 days, and Tcl = 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 casesa
Clinical case number
Tumor doubling time T d (0)
Cell clearance time T cl
Retardation factor θ
Irradiation day for GKSRS
Tumor volume at GKSRS
Tumor doubling time at GKSRS
Volume ratio, R40 b
% Cell survival fraction
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.
The tumor volume doubling time at the time of GKSRS, Td(GKSRS), ranged from 7.9 days to 32.6 days.
The cell clearance time is between 10 and 40 days.
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.
P-value for predictors of response R40 a
P-value: 2nd order
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, Td(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 . 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 well-known 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 . 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 non-dimensional 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 petri-dish 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, Tm, 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 non-dividing cells, i.e., sublethal-cell 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 . 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 self-renewal 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 MRI-based 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 two-component (or two-compartment) model that was developed in the 1970s . 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 non-dividing 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 head-and-neck cancer cases [15–17]. We also adopted this cell clearance model in this study.
Recently, the two-component model was extended to a three-component model that includes a reversible state named as “cell-cycle arrest” state and denoted by C  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 in-vitro 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 ; 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 . 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 non-uniform nature of the dose distributions. This warrants further investigation.
Delineation of tumor volume
In this study, we assumed that the volume visible as a contrast-enhanced area on T1-weighted 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 . 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 non-dividing 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.
Our model contains eight biological parameters, i.e., α, α/β, Td (tumor doubling time), Tcl (cell clearance time), θ (vascular growth retardation factor), τrad (effective time of radiation), Tcc (cell cycle time), and Tm (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, R40, 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 R40 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 patient-specific 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.
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 radiation-damaged 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 post-GKSRS tumor volume change and the α and θ-values.
Further refinement of the model which includes the radiation-induced 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.
All calculations were executed on a PC using MATLAB. The code is available from the first author.
contrast enhanced T1-weighted magnetic resonance imaging
cancer stem cells
Gamma Knife stereotactic radiosurgery
magnetic resonance imaging
mean square of differences
non-small cell lung
ordinary differential equation
renal cell carcinoma
The part of this work was previously presented at the American Association of Physicists in Medicine (AAPM) annual meetings in the poster format, SU-ET-08 and SU-E-T-309, 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.
- 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 J-P, 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/s00285-008-0215-x.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. Single-cell-based computer simulation of the oxygen-dependent tumour response to irradiation. Phys Med Biol. 2007;52(16):4775–89. doi:10.1088/0031-9155/52/16/005.View ArticlePubMedGoogle Scholar
- Titz B, Jeraj R. An imaging-based 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/s00285-008-0219-6.PubMed CentralView ArticlePubMedGoogle Scholar
- Perez-Garcia VM, Bogdanska M, Martinez-Gonzalez A, Belmonte-Beitia J, Schucht P, Perez-Romasanta LA. Delay effects in the response of low-grade gliomas to radiotherapy: a mathematical model and its therapeutical implications. Math Med Biol. 2014. doi:10.1093/imammb/dqu009
- Nawrocki S, Zubik-Kowal 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/0008-5472.CAN-09-2501.PubMed CentralView ArticlePubMedGoogle Scholar
- Chvetsov AV. Tumor response parameters for head and neck cancer derived from tumor-volume 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. Tumor-volume simulation during radiotherapy for head-and-neck cancer using a four-level 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 X-irradiated rat embryo cells transfected with c-myc or c-Ha-ras. 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 x-rays. 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 x-irradiation. Eur J Cancer. 1969;5(2):173–89.View ArticlePubMedGoogle Scholar
- Thompson LH, Suit HD. Proliferation kinetics of x-irradiated mouse L cells studied WITH TIME-lapse 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 low-grade glioma treated with chemotherapy or radiotherapy. Clin Cancer Res. 2012;18(18):5071–80. doi:10.1158/1078-0432.CCR-12-0084.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 x-rays 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 H-Y, 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/blood-2008-08-077941.View ArticlePubMedGoogle Scholar
- Li L, Bhatia R. Stem cell quiescence. Clin Cancer Res. 2011;17(15):4936–41. doi:10.1158/1078-0432.CCR-10-1499.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 PDGF-driven 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
- Garcia-Barros M, Paris F, Cordon-Cardo C, Lyden D, Rafii S, Haimovitz-Friedman 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