A mathematical model of tumor growth and its response to single irradiation

Abstract

Background

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.

Methods

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).

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 post-treatment 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 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 [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 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 [13] 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. [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 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.

Mathematical model

The proposed tumor-volume 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 two-components, in the tumor, dividing cells, and non-dividing 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 tumor-growth rate λ(t) times the cell proliferation probability p(D), or can transit to the non-dividing state with a rate g(D). The non-dividing cells are eventually cleared from the original location with the cell clearance rate η_cl.

Fig. 1

Diagram of the proposed 2-component model

This model can be represented by a set of three equations of three variables: the volume of proliferating tumor, V_T(t), the volume composed of non-dividing 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:

$$ \frac{d{V}_T}{dt}=\lambda (t){\displaystyle {\boldsymbol{V}}_T} $$

(1)

$$ \frac{d\lambda }{dt}=-\theta \lambda (0)\lambda $$

(2)

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 radiation-effect time τ_rad, i.e., in t_R ≤ t < t_R + τ_rad:

$$ \frac{d{V}_T}{dt}=\lambda (t)p(D){\displaystyle {\boldsymbol{V}}_T}-g(D){\displaystyle {\boldsymbol{V}}_T} $$

(3)

$$ \frac{d{\displaystyle {V}_{ND}}}{dt}=g(D){\displaystyle {V}_T}-{\displaystyle {\eta}_{cl}}{\displaystyle {V}_{ND}} $$

(4)

$$ \frac{d\lambda }{dt}=-\theta \lambda (0)\lambda $$

(5)

For t > t_R + τ_rad, V_T, V_ND, and λ are the solutions of the following ODEs:

$$ \frac{d{V}_T}{dt}=\lambda (t){\displaystyle {\boldsymbol{V}}_T} $$

(6)

$$ \frac{d{\displaystyle {\boldsymbol{V}}_{ND}}}{dt}=-{\displaystyle {\boldsymbol{\eta}}_{cl}}{\displaystyle {\boldsymbol{V}}_{ND}} $$

(7)

$$ \frac{d\lambda }{dt}=-\theta \lambda (0)\lambda $$

(8)

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 non-dividing, g(D), are given as a function of the characteristic time T*, the observation time T_m, and a dose-dependent constat χ(D). Furthermore, for the LQ model, χ(D) is a function of α and α/β and it is given by [29]

$$ \chi (D)=\alpha D\left(1+\frac{D}{\alpha /\beta}\right) $$

(9)

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 non-dividing cells is V _ND . The non-dividing cells are removed from the original site. These processes are represented by the following two ordinary differential equations (ODEs).

$$ \frac{d{V}_T}{dt}=f\left({\displaystyle {V}_T},{\displaystyle {V}_v},D\right){\displaystyle {V}_T}-g(D){\displaystyle {V}_T} $$

(10)

$$ \frac{d{\displaystyle {V}_{ND}}}{dt}=g(D){\displaystyle {V}_T}-{\displaystyle {\eta}_{cl}}{\displaystyle {V}_{ND}} $$

(11)

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 [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³. 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

$$ f\left({\displaystyle {V}_T},{\displaystyle {V}_v},D\right)=\lambda \left({\displaystyle {V}_T},{\displaystyle {V}_v}\right)p(D) $$

(12)

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 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:

$$ \lambda =a\frac{{\displaystyle {V}_v}}{{\displaystyle {V}_T}} $$

(13)

By taking the time derivative of Equation (13), one obtains

$$ \frac{d\lambda }{dt}=\left\{\frac{1}{{\displaystyle {V}_v}}\frac{d{\displaystyle {V}_v}}{dt}-\frac{1}{{\displaystyle {V}_T}}\frac{d{\displaystyle {V}_T}}{dt}\right\}\lambda $$

(14)

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

$$ \frac{d\lambda }{dt}=\left\{{\displaystyle {\lambda}_v}(0)-\lambda (0)\right\}\lambda $$

(15)

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:

$$ {\displaystyle {\lambda}_v}(t)=\theta \lambda (t) $$

(16)

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

The standard model for radiation-induced 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 linear-quadratic (LQ) equation for the cell survival fraction S, then, the volumes of the dividing cells and the non-dividing cells after a single pulse of irradiation at time t ₊ jump from those at pre-irradiation time t ₋:

$$ {\displaystyle {V}_T}\left({\displaystyle {t}_{+}}\right)=S(D){\displaystyle {V}_T}\left({\displaystyle {t}_{-}}\right) $$

(17)

$$ {\displaystyle {V}_{ND}}\left({\displaystyle {t}_{+}}\right)=\left\{1-S(D)\right\}{\displaystyle {V}_T}\left({\displaystyle {t}_{-}}\right) $$

(18)

where S(D) = e ^−χ(D) and χ(D) = αD + βD ². With this model, the function f and g do not change before and after irradiation and those are given by

$$ f=\lambda $$

(19)

$$ g=0 $$

(20)

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 non-dividing 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 (=1-p) 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 non-dividing cells, N _ND :

$$ \frac{d{\displaystyle {N}_D}}{dt}=\frac{2p-q}{T^{*}}{\displaystyle {N}_D} $$

(21)

$$ \frac{d{N}_{ND}}{dt}=\frac{q}{T^{*}}{N}_D $$

(22)

Here T ^* is a characteristic time, which can be set equal to T _cc . The initial conditions are

$$ {\displaystyle {N}_D}(0)={\displaystyle {N}_0} $$

(23)

$$ {\displaystyle {N}_{ND}}(0)=0 $$

(24)

The solutions of Equations (21) and (22) can be easily found as follows:

$$ {\displaystyle {N}_D}(t)={\displaystyle {N}_0} \exp \left(\frac{2p-q}{T^{*}}t\right) $$

(25)

$$ {\displaystyle {N}_{ND}}(t)=\frac{q}{2p-q}{\displaystyle {N}_0}\left\{ \exp \left(\frac{2p-q}{T^{*}}t\right)-1\right\} $$

(26)

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:

$$ \exp \left(\frac{2p-q}{T^{*}}{\displaystyle {T}_m}\right) \exp \left(-\frac{2}{T^{*}}{\displaystyle {T}_m}\right)=S(D)= \exp \left(-\chi (D)\right) $$

(27)

By solving Equation (27) for p, we obtain

$$ p(D)=1-\frac{T^{*}}{3T{}_m{}}\chi (D) $$

(28)

The function g(D) is equal to q/T _m :

$$ g(D)=\frac{\chi (D)}{3{T}_m} $$

(29)

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 radiation-effect 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

$$ {\displaystyle {\int}_{{\displaystyle {t}_{-}}}^{{\displaystyle {t}_{-}}+{\tau}_{rad}}g(D)}{\displaystyle {N}_D}(t)dt\approx \left\{1-S(D)\right\}{\displaystyle {N}_D}\left({\displaystyle {t}_{-}}\right) $$

(30)

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 non-dividing 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 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.

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₁,…, 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:

$$ {\displaystyle \underset{p\in \varGamma }{ \min }}{\displaystyle {\sum}_{i=1}^N{\displaystyle {\left\{{\displaystyle {y}_i}-{\widehat{y}}_i\left(\overrightarrow{p}\right)\right\}}^2}} $$

(31)

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 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

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

$$ MSD = \frac{1}{N}{\displaystyle \sum_{n=1}^N}{\left(\frac{y_n-{Y}_n}{Y_n}\right)}^2 $$

(32)

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:

$$ AIC=N\bullet ln(MSD)+2k $$

(33)

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 X-ray or neutron beam. The data are reproduced as Figure 20.2A in the radiobiology textbook [29]. 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 [37]. 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 [38]. 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.

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., non-small 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].

Table 1 Treatment-related parameters for clinical Gamma Knife stereotactic surgery cases

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₄₀, 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₄₀ 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 non-linear regression of the second-order polynomial function between R₄₀ 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.

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.

Table 2 Model parameters for rat rhabdomyosarcoma experiment data: the proposed model^a

Fig. 2

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 proposed 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⁻¹, θ = 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.

Fig. 3

Tumor volume as a function of time after the first diagnostic MRI scan. The time is in days. The circles indicate the tumor volume measured on Gd contrast enhanced MRI of GK patients. Solid lines are the volume variation predicted by the model. The vertical arrows indicate the time of GKSRS. The vertical arrows indicate the time of GKSRS. The error bars of the measured volumes were obtained by calculating the volume of an equivalent sphere with ±1 mm radius of the original volume. a case one, b case two, c case three, d case four, and e case five

Table 3 Model parameters for radiation response: clinical GKSRS cases^a

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 P-values of the correlation between R₄₀ 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₄₀ and α or θ.

Table 4 P-value for predictors of response R₄₀ ^a

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 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 [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 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, 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 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 [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 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 [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 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 [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 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 [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⁻¹), which are on the smaller side of the commonly accepted range between 0.1 and 0.3 Gy⁻¹. 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 [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 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.

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₄₀, 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₄₀ 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.

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 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.

Methods

Computational calculations

All calculations were executed on a PC using MATLAB. The code is available from the first author.