 Research
 Open Access
 Published:
A mathematical model of tumor growth and its response to single irradiation
Theoretical Biology and Medical Modelling volume 13, Article number: 6 (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.
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. 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
In this section, we summarized the proposed model used in this study. A detailed description of the derivation of the equations and biologicl assumptions are given in the following sections. Two cell types, or twocomponents, in the tumor, dividing cells, and nondividing cells, were included in the model. Figure 1 illustrates that under the influence of radiation dose D, the dividing cells can continue dividing with the rate of the tumorgrowth rate λ(t) times the cell proliferation probability p(D), or can transit to the nondividing state with a rate g(D). The nondividing cells are eventually cleared from the original location with the cell clearance rate η_{cl}.
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}.
Initially, the tumor volume grows, according to the following ODEs:
After a single instantaneous irradiation of dose D given to the tumor at t_{R}, the following ODEs are solved for a time period of the active radiationeffect time τ_{rad}, i.e., in t_{R} ≤ t < t_{R} + τ_{rad}:
For t > t_{R} + τ_{rad}, V_{T}, V_{ND}, and λ are the solutions of the following ODEs:
Note that three functions, V_{T}(t), V_{ND}(t), and λ(t), are continuous at t_{R} and t_{R} + τ_{rad}. The probability p(D) and the transition rate of tumor cells from dividing to nondividing, g(D), are given as a function of the characteristic time T*, the observation time T_{m}, and a dosedependent constat χ(D). Furthermore, for the LQ model, χ(D) is a function of α and α/β and it is given by [29]
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
By assuming that the volume of tumor is proportional to the number of cells in the tumor [30], we considered the volume as the basic parameter instead of the tumor mass. We denoted the total volume of malignant cells by V _{ T }. A malignant tumor cell divides with a specific cell cycle time, T _{ cc }. When these cells are exposed to radiation, some of them stop dividing and die. Some fraction of the irradiated cells survives and keeps dividing. The volume of the nondividing cells is V _{ ND }. The nondividing cells are removed from the original site. These processes are represented by the following two ordinary differential equations (ODEs).
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
Tumors require new vasculature to provide additional nutrient to grow beyond a certain size, i.e., about 1 mm^{3}. Angiogenesis, the process by which vascular structure develops inside the tumor, is the main mechanism to keep the tumor supplied with sufficient nutrients while it is growing. We represented the amount of the nutrient supplied by the volume of blood or the vascular structure in the tumor, V _{ v }. Then, we assumed that the function f in Equation (10) is made of two factors and it is given by
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 growth rate of the tumor may depend on the blood volume. Hence, we assumed that the growth rate λ is proportional to the ratio of the vascular and tumor volumes:
By taking the time derivative of Equation (13), one obtains
This indicates that the rate of change in λ is the difference between the specific volume growth rates of the vasculature and the tumor. For simplicity, we replaced these growth rates with the initial values. Then, we have
We assumed that the volume of the vascular structure increases with the growth rate of λ_{v} as the tumor grows, but the vascular volume grows slower than the tumor:
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
The standard model for radiationinduced cell death after a single dose assumes that a portion of the tumor cells dies and the rest still maintains its proliferation capacity [12, 15, 16, 20, 32]. By applying the linearquadratic (LQ) equation for the cell survival fraction S, then, the volumes of the dividing cells and the nondividing cells after a single pulse of irradiation at time t _{+} jump from those at preirradiation time t _{−}:
where S(D) = e ^{−χ(D)} and χ(D) = αD + βD ^{2}. With this model, the function f and g do not change before and after irradiation and those are given by
Furthermore, V _{ ND }(t) = 0 for t ≤ t _{−} for single fraction treatment.
Some studies have shown that cancer cells do not lose their mitotic capability right after exposure to radiation, but those may continue dividing [21–24]. Hence, the volumes of the dividing cells and nondividing cells gradually change in time. The following radiation response model was formulated by considering this phenomenon. Individual cell divides with probability p and does not divide with probability q (=1p) in the cell cycle time T _{ cc }. Then, it is easy to derive the following rate equations for the number of dividing cells, N _{ D }, and the number of nondividing cells, N _{ ND }:
Here T ^{*} is a characteristic time, which can be set equal to T _{ cc }. The initial conditions are
The solutions of Equations (21) and (22) can be easily found as follows:
We applied these expressions to a radiobiological experiment [33], by which one measures the cell survival fraction after a single dose of radiation, D. The measured survival fraction S(D) is the ratio of the number of dividing (or clonogenic) cells at a specific observation time, or colony counting time, T _{ m }, when the cells are irradiated and the number of dividing cells at T _{ m } when the cells are not irradiated. The data for the surviving fraction is fitted well by the LQ equation [23, 29, 34]. Hence, we obtain the following relationship:
By solving Equation (27) for p, we obtain
The function g(D) is equal to q/T _{ m }:
Note that Equations (28) and (29) indicate the probabilities p(0) = 1 and q(0) = 0 when D = 0 as expected. For this study, we assumed that the dose is given by a single fraction in duration negligible relative to the time scale in which we monitor the volume variation. Equation (28) indicates that the mitotic probability of the cell decreases from unity to a value smaller than unity because of radiation. At the same time, the cells die at the rate of g(D). It was assumed that the effect of the radiation on the cell killing probability lasts for only a certain length, active radiationeffect time denoted as τ _{ rad }. In other words, after the period τ _{ rad }, the probabilities p and q are back to 1.0 and 0.0, respectively. It is noteworthy that the total number of proliferating cells dying after irradiation is approximately related to the survival fraction S(D) by
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.
The parameter estimation was also automatically performed. To fit the mathematical model to the measured tumor volume, we combined the least squares approach with the simulated annealing algorithm [35]. In particular, we needed to compare the N observed tumor volumes denoted by (y_{1},…, y_{N}) and the model predicted volumes at the corresponding times. For a set of parameters \( \overrightarrow{p} = \left(\alpha, \uptheta,\ {\mathrm{T}}_{\mathrm{d}},\ {\mathrm{T}}_{\mathrm{cl}}\right) \), we denoted the model predicted volumes at the corresponding times by \( \left({\widehat{y}}_1\left(\overrightarrow{p}\right),\dots, {\widehat{y}}_N\left(\overrightarrow{p}\right)\right) \). Then our method for finding the best value of \( \overrightarrow{p} \) was to solve the minimization problem:
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
To evaluate the model performance for data fitting of the experimental data, we calculated the mean square of differences (MSD) between the experimental data \( \left\{{Y}_n\right\} \) and the volumes estimated by the models at the corresponding time \( \left\{{y}_n\right\} \). MSD was defined by calculating a relative error at each measurement point. For N measurement points of data, MSD was given by
In addition to MSD, the quality of models can be compared by calculating the Akaike Information Criterion (AIC) [36]. We used the following definition of AIC:
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.
Five clinical cases analyzed for this study are summarized in Table 1. The local control of the tumor after GKSRS are categorized into four groups: complete response (CR), partial response (PR), stable disease (SD), and progressive disease (PD) [39]. If the size of the tumor does not change more than 20 % after the treatment, it is designated by SD. CR means the disappearance of the tumor. The volume change between CR and SD is considered PR. PD indicates the treatment failure. Using this definition, the local control status of five cases 1 to five was assigned to CR, PR, PR, SD, and PD, respectively. The complete response which was seen in case one could be mainly due to the large mean dose (38.3 Gy) given to the tumor in addition to other factors such as the primary tumor type, i.e., nonsmall cell lung cancer. The failure of the local control of case five could be explained by the low maximum dose and the type of the primary tumor. In fact, there is a clinical evidence showing that 20 Gy marginal dose is not sufficient for melanoma metastasis for local control [40].
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
The model parameters of the rat rhabdomyosarcoma experiments are given in Table 2. The initial tumor doubling time, T_{d}(0), the α/β value, and the cell clearance time, T_{cl}, were set to the same values for the five cases. The quality of the parameter fitting can be seen in Fig. 2. The experimental data are shown as points in the figure, whereas the volume changes estimated by the model are plotted as solid lines for various dose levels, i.e., 0 (or control), 1000, 2000, 3000, and 4000 rads. The excellent data fitting was achieved by adjusting only the α value and the retardation factor θ for five different radiation doses. To quantify the magnitude of radiation effect, the survival fraction, S, of the dividing cells at the end of the active radiation period τ _{ rad } was the lower and the values are shown in the lowest row in Table 2. It is noted that the plotting of S as the function of dose generates the typical curve obtained by the LQ model.
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
The clinically observed changes of tumor volumes are presented as open circles in Fig. 3a–e. The solid lines in the figure indicate the temporal volume change generated by the solutions of the mathematical model developed for this study. Four parameters, α, T _{ d } (0), T _{ cl } , and θ were optimized for the best fit. Table 3 shows the model parameters used to obtain these curves. Note that the radiobiological parameter α/β = 10 Gy, τ _{ rad } = 8 days, and the cell cycle time T_{cc} = 1 day for all cases. We assumed that the number of cell cycles in which the survival fraction was measured to be ten, or T_{m} = 10 days. The dose used for the simulation was the mean dose prescribed to the target. Due to the technical limitation, the measured volume suffers from a large uncertainty. To assess the uncertainty of the measured tumor volume, hence, we calculated the volume of an equivalent sphere with 0.1 cm expansion or reduction of the radius of the sphere having the same volume as the original volume of the irregularly shaped tumor. Those values are indicated as error bars in Fig. 3.
From the model parameters shown in Tables 1 and 3 and the tumor growth patterns seen in Fig. 3, we can make several observations:

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.
The results of the linear regression analyses were shown in Table 4. The Pvalues of the correlation between R_{40} and five variables, T_{d}(GKSRS), α, θ, V_{tumor}(GKSRS), and D_{mean}, were 0.302, 0.026, 0.007, 0.324, and 0.817, respectively. The result indicates a significant correlation between R_{40} and α or θ.
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
References
 1.
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.
 2.
Barillot E, Calzone L, Hupe P, Vert JP, Zinovyev A. Computational systems biology of cancer. Boca Raton: CRC Press; 2013.
 3.
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.
 4.
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.
 5.
Deisboeck TS, Stamatakos GS. Editors. Multiscale cancer modeling. Chapman and hall/CRC mathematical and computational biology (book 34). Boca Rayton: CRC Press; 2010.
 6.
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.
 7.
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.
 8.
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.
 9.
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.
 10.
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.
 11.
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.
 12.
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.
 13.
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
 14.
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
 15.
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.
 16.
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.
 17.
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.
 18.
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.
 19.
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.
 20.
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.
 21.
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.
 22.
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.
 23.
Joiner M, van der Kogel A, editors. Basic clinical radiobiology. 4th ed. London: Hodder Arnold; 2009.
 24.
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.
 25.
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.
 26.
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.
 27.
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.
 28.
Tannock I, Howes A. The response of viable tumor cords to a single dose of radiation. Radiat Res. 1973;55(3):477–86.
 29.
Hall EJ, Giaccia AJ. Radiobiology for the radiologist. 7th ed. Philadelphia: Lippincott Williams and Wilkins; 2011.
 30.
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.
 31.
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.
 32.
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.
 33.
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.
 34.
Dale RG, Jones B, editors. Radiobiological modelling in radiation oncology. Oxfordshire: The British Institute of Radiology; 2007.
 35.
Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical recipes: the Art of scientific computing. 3rd ed. New York: Cambridge University Pres; 2007.
 36.
Brockwell PJ, Davis RA. Time series: theory and methods. 2nd ed. New York: Springer; 1991.
 37.
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.
 38.
Greenwood J. Mechanisms of blood–brain barrier breakdown. Neuroradiology. 1991;33(2):95–100. doi:10.1007/BF00588242.
 39.
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.
 40.
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.
 41.
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.
 42.
Dick JE. Stem cell concepts renew cancer research. Blood. 2008;112(13):4793–807. doi:10.1182/blood200808077941.
 43.
Li L, Bhatia R. Stem cell quiescence. Clin Cancer Res. 2011;17(15):4936–41. doi:10.1158/10780432.CCR101499.
 44.
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.
 45.
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.
 46.
Roberts TP, Chuang N, Roberts HC. Neuroimaging: do we really need new contrast agents for MRI? Eur J Radiol. 2000;34(3):166–78.
 47.
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.
 48.
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.
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).
Author information
Affiliations
Corresponding author
Additional information
Competing interests
All authors declare that they have no competing interests.
Authors’ contributions
YW conceived the ideas, wrote the MATLAB codes, carried out the analyses, and prepared the initial manuscript. ED collected the clinical data used in this study. KL wrote the MATLAB code and did parameter estimation analysis using the simulated annealing technique. SH proposed using the rat experimental data for the study. All reviewed the manuscript. All authors read and approved the final manuscript.
Additional files
Additional file 1:
Mathematical details. (DOCX 27 kb)
Additional file 2:
Comparison of the proposed model with a traditional model. Table S1. Model parameters for rat rhabdomyosarcoma experiment data: the old model^{(*)}. Figure S1. Temporal change of the rat rhabdomyosarcoma tumor volume before and after irradiation. The discrete points indicate the experimental data, whereas the solid lines show the tumor volume change produced by the old model. (DOCX 174 kb)
Rights and permissions
Open Access This 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.
About this article
Cite this article
Watanabe, Y., Dahlman, E.L., Leder, K.Z. et al. A mathematical model of tumor growth and its response to single irradiation. Theor Biol Med Model 13, 6 (2016). https://doi.org/10.1186/s1297601600327
Received:
Accepted:
Published:
Keywords
 Mathematical model
 Tumor volume
 Gamma knife
 Radiation effect