Kinetic modeling of tumor regression incorporating the concept of cancer stem-like cells for patients with locally advanced lung cancer.

BACKGROUND
Personalized medicine for patients receiving radiation therapy remains an elusive goal due, in part, to the limits in our understanding of the underlying mechanisms governing tumor response to radiation. The purpose of this study was to develop a kinetic model, in the context of locally advanced lung cancer, connecting cancer cell subpopulations with tumor volumes measured during the course of radiation treatment for understanding treatment outcome for individual patients.


METHODS
The kinetic model consists of three cell compartments: cancer stem-like cells (CSCs), non-stem tumor cells (TCs) and dead cells (DCs). A set of ordinary differential equations were developed to describe the time evolution of each compartment, and the analytic solution of these equations was iterated to be aligned with the day-to-day tumor volume changes during the course of radiation treatment. A least squares fitting method was used to estimate the parameters of the model that include the proportion of CSCs and their radio-sensitivities. This model was applied to five patients with stage III lung cancer, and tumor volumes were measured from 33 cone-beam computed tomography (CBCT) images for each of these patients. The analytical solution of these differential equations was compared with numerically simulated results.


RESULTS
For the five patients with late stage lung cancer, the derived proportions of CSCs are 0.3 on average, the average probability of the symmetry division is 0.057 and the average surviving fractions of CSCs is 0.967, respectively. The derived parameters are comparable to the results from literature and our experiments. The preliminary results suggest that the CSC self-renewal rate is relatively small, compared to the proportion of CSCs for locally advanced lung cancers.


CONCLUSIONS
A novel mathematical model has been developed to connect the population of cancer stem-like cells with tumor volumes measured from a sequence of CBCT images. This model may help improve our understanding of tumor response to radiation therapy, and is valuable for development of new treatment regimens for patients with locally advanced lung cancer.


Backgrounds
Radiation therapy (RT) is the treatment of choice for many lung cancer patients. The aim of an RT plan is to deliver high radiation dose to tumor volume to eliminate cancer cells without causing severe damage to surrounding healthy tissue. The efficacy of a planned treatment, often quantified by tumor control probability (TCP), depends on several patient-specific factors such as the surviving fraction of tumor cells and treatment regimens used [1,2]. Many models and methods have been proposed to calculate the control probability [3]. In general, TCP is defined as an exponential function of the number of surviving cells after a given treatment [4]. The number of the surviving cells depends on the radiation dose and the radio-sensitivity of each cell subpopulation within the tumor. Much effort has been devoted to the development and validation of treatment response models such as, for example, the linear quadratic and the modified linear quadratic models [5,6]. However, without calibrating radio-sensitivity parameters such as the α/β ratio for individual patients [7,8], inter-patient differences cannot be taken into account in these models [9][10][11].
The radio-sensitivity of tumor cells could be heterogeneous among its subpopulations [12]. Bhaumik and Jain divided cells into the proliferating and non-proliferating compartments, and investigated the tumor growth kinetics on the basis of interactions between the two compartments [13]. Rockne et al. extended the reactiondiffusion model of tumor proliferation and invasion to a general form to include tumor response to radiation therapy [14,15], and Ribba et al. proposed a multiscale model of cancer to investigate tumor response to radiation [16]. Zhong et al. used a pair of ordinary differential equations to capture the essential features of tumor dynamics, modeling radiation-induced kinetics in the number of tumor cells [17,18]. Zaider showed that the amount of radiation dose necessary for obtaining a 90% TCP can be dominated by radio-resistant cells which could be as small as a millionth of the total number of tumor cells [11]. The challenge to all these approaches has been the lack of a measure of the number of these cells and their radio-sensitivities [9,19,20].
In the context of cancer stem-like cells (CSCs) which, by definition, are self-renewing with unlimited proliferative potential and a capacity to generate differentiated cells (i.e. cells that evolve to express proteins that confer advantages for growth, metastasis, invasion, and resistance to treatment) [21], Hillen et al. developed a mathematical model of a heterogeneous population of CSCs and non-stem tumor cells (TCs) to investigate movement between these components [19]. With this model, Yu et al. investigated the impact of different radio-sensitivities between stem-like and non-stem cancer cells, and demonstrated the importance of CSCs for tumor local control [22]. In these studies, the CSC self-renewal rate was assumed to be 1% [19,22,23]. These models have helped us to get a better understanding of the mechanism of heterogeneous tumor response to radiation treatment. However, the parameters used in these models are associated with large uncertainties [24]. For example, the population of CSCs has been reported to range from 0.1 to 82.5% [25][26][27]. It has been difficult to identify CSCs and measure their subpopulations and radio-sensitivity for individual patients. The lack of an efficient method to measure these radiobiological parameters has limited the clinical use of these models for treatment planning and monitoring response in individual patients.
With advances in radiation therapy techniques, it has become possible to measure changes in tumor volume, a major index of treatment response for solid tumors, using daily anatomical images acquired during the course of radiation treatment. In this study, we have developed a set of differential equations that includes the compartments of cancer stem-like and non-stem cells, and dead cells so that changes in each of these cell compartments affecting tumor volume during treatment can be evaluated. The analytical solution of these equations, after validating using numerical simulations, was applied to treatment regimens to form an iterative model for individual patients. Based on the tumor volumes measured from daily CBCT images, patient-specific parameters including the proportion of CSCs and their radio-sensitivities were estimated, and stabilities of these parameters were assessed. In the next sections the theoretical model is developed, and applied to a sequence of CBCT images to evaluate radiobiological parameters. The derived parameters are compared with those reported in the literature.

Materials and methods
The observations of radiosensitivity of normal tissue stem cells by Bergonie and Tribondeau at the turn of the last century, confirmed by 100 years of experience with radiation, indicate that normal tissue stem cells are typically more sensitive to radiation than differentiated cells [28]. However, many factors such as cell quiescence, DNA repair ability, overexpression of anti-apoptotic proteins, detoxifying enzymes such as glutathione and a hypoxic niche microenvironment render CSCs to be radiation resistant and responsible for tumor growth and relapse [29][30][31]. CSCs shared some characteristics in common with normal tissue stem cells [32]. For example, CSCs were assumed to exhibit unlimited proliferation and a low baseline apoptosis [33]. From CSCs to differentiated cells, there are multiple stages of differentiations. Due to the limited availability of clinical and experimental data, it is impossible to model cellular interactions at all these stages. In this study we divide cells in a tumor volume observable in CBCT images into three compartments: cancer stem-like cells, non-stem tumor cells (TC) and dead cells (DC).

A kinetic model of tumor response to radiation treatment
To model the transient status of hierarchical cancer cells, each cell subpopulation is described as a compartment, and different compartments are connected by transitional events with fixed rates. For example, a cancer stem-like cell is assumed to divide at a constant rate, to either two CSCs (symmetric division) or one CSC and one TC (asymmetric division). Let u(t) and v(t) denote the numbers of CSCs and TCs per unit volume, and w(t) the number of DCs before being moved out of the tumor volume. Let m u and m v denote the division rates of CSCs and TCs. Let δ represent the probability of the symmetry division which is equivalent to the rate of the CSC self-renewal, and let ϒ u and ϒ v denote the rates of the radiation-induced cell death for CSCs and TCs, m a the rate of the programmed cell death for TCs without irradiation, and ϒ w the clearance rate of DCs. The interaction of the three compartments (CSC, TC and DC) can be mathematically represented by To apply this set of differential equations to daily treatment regimens which specifies when the treatment will be delivered and how much radiation dose will be used during the course of radiation therapy, these equations have to be solved explicitly. Let SF u and SF v denote the surviving fraction of cells in u(t) and v(t) after being irradiated by one daily fraction, 2 Gy. ϒ u and ϒ v can be represented as ϒ u = − ln SF u and ϒ v = − ln SF v . The parameters m u , m v are related to the tumor doubling times T u and T v by m u = ln 2/T u and m v = ln 2/T v , and the clearance rate of DCs is ϒ w = ln 2/T h , and m a = ln 2/T a where T h is the half time for the dead cell recycling and T a is the half time for the programmed cell death. Following the mathematical reasoning illustrated in [17], the differential equations can be explicitly solved as Let t = 1, then the resultant solution can be applied to each RT fraction or non-RT days: for non-RT days SF u and SF v will be set to 1; otherwise, these parameters will be optimized. The resultant kinetic model is executed on each day, and the computed numbers of the cells in each compartment at day i + 1 can be written as: Since CBCT imaging was used in patient treatment, radiation impact on cell survival can be assessed on a daily basis. By iterating Eq. (3) for each day, the number of the survived cells in each compartment can be calculated, and consequently the number of the total survived cells can be compared with the CBCT-measured results.

The numerical solution of the kinetic model
With any given parameter, Eq. (1) can be numerically solved. The numerical solution can be used to validate those generated by the iterative functions in Eq. (3). Specifically, an Euler method was applied to Eq. (1) with the procedure sketched below: Let X(t) = (u(t), v(t), w(t)), and let _ X t ¼ f ðt; XÞ with f(t, X) defined by the differential Eq. (1). Let the initial condition satisfy X(0) = X 0 = (u 0 , v 0 , 0). For k = 0, . . ., n-1, let a = k, b = k + 1, and h = (ba)/m. Then for i = 0, …, m, The vector X m will be an approximation to the solution X(k + 1) of the differential equations, for each k = 0,. .., n-1. By adding all the components in X m together, the number of the total tumor cells can be calculated and compared to the results derived from the analytical iterative model in Eq. (3).
Equation (4) can be used to validate the analytical solution (3), but is not suitable to be used for parameter optimization. Note that the numerical method costs more computation time than using the explicit, analytic solution, and also the limited precision of the numerical solution may compound the uncertainties of the parameter derivation. After verification of the consistency of the two models, only the analytical model in (3) was used to derive radio-biological parameters. The simulation of treatment regimens and its related parameter optimization were programmed using C++ on a Linux workstation.

Tumor volume measurement from CBCT images for lung cancer patients
With the approval from the Institutional Review Board (IRB) of Henry Ford Health System (HFHS), five lung cancer patients treated at HFHS were investigated in this study. Among the five patients, 2 are with adenocarcinoma and 3 with squamous cell cancer (see Table 1). All of them were at stage III, and treated 33 fractions, 2 Gy/ fraction using 6 MV IMRT. Daily CBCT images were acquired to help setup patients and monitor tumor regression. Gross tumor volumes (GTVs) were contoured by physicians on each of these CBCT images, and the contoured volumes were calculated using the Eclipse planning system (Varian Medical Systems, Palo Alto, CA). The dates for RT and CBCT imaging, and the tumor volume measured on each of these images were written into an input file for the kinetic model. Most of the patients have their treatment course completed within 47 days, except patient A2 who had no treatment for one week, resulting in the RT course extended to 55 days.
A least squares fitting method to optimize the model For the five patients investigated in this study, the tumor volumes delineated in 33 CBCT images were normalized to that in the first CBCT to obtain the rate of tumor regression for each of these RT days. The coefficients δ, p, T h , and SF u in the model were estimated by comparing the modeled tumor volumes with those measured from CBCT images using a least squares fitting method. With the assumption that tumor volume is proportional to the number of its clonogens with a linear constant coefficient [34], an objective function can be defined by where M i is the tumor volume measured from the i-th CBCT scan, and G i = U i + V i + W i is its corresponding cell number predicated by the kinetic model (3). At the beginning of treatment, it is assumed that W 0 = 0 and the proportion of CSCs is denoted by p = U 0 /(U 0 + V 0 ). Let T u , T v , and T a share the same doubling time as suggested in [22][23][24], and the surviving fraction SF T of the total cells be written as The two terms in the right-hand side of Eq. ( (6)) represent the survival fractions of stem and non-stem cells, respectively. With the parameters T u and SF T measured for different types of cancer cells, the other parameters were derived by minimizing the function R in Eq. ((5)) using the least squares fitting method.
It has been reported that the cell surviving fraction under 2 Gy is 0.90 and the cell division time T div is 60 h for squamous cell cancer [35,36]. For adenocarcinoma (A549), it has been reported that the surviving fraction is 0.62 and cell division time is 22 h [37,38]. However, T div is not appropriate to be used directly as the tumor doubling time in Eq. (3) due to the existence of quiescent cells. On the other hand, the time (T vol ) required to double tumor volume was reported to be 166.3 days on average over 237 patients, with 221.6 days for adenocarcinoma, 115.2 days for squamous cell carcinoma [39]. Note that the measured time (T vol ) is based on the resultant volumes that already included the factor of cell loss, and therefore, does not represent the cell division time required in Eq. (3). Different from the above two rates of tumor growth, the potential doubling time (T pot ), calculated as the ratio of the time for DNA synthesis and the proportion of cells synthesizing DNA, takes into account the number of quiescent cells, and is a more objective index estimating proliferation occurring during a course of radiation treatment [40]. Typical values for T pot have been found to be 6.2 days for squamous and 7.1 days for adenocarcinoma cells [22,41]. In the next section, we will compare the potential doubling time with the cell division and volume doubling times, and evaluate the impact of tumor growth uncertainties on the optimized parameters. We will also measure the cell survival fraction and compare the measured results with those reported in literature for squamous cell carcinoma.

Results of the kinetic modeling
Based on the potential doubling times reported for adenocarcinoma and squamous cell cancer, the least squares fitting method was applied to the CBCT image sequences acquired from the five patients to minimize the objective function R. The resultant residues were found to be less than 8.2%. The optimized parameters were listed in Table 2 With these optimized parameters, the number of tumor cells in each compartment was calculated using the iterative function in Eq. (3). The total number of the tumor cells computed for each day was normalized to that on the first RT day to get compared with the relative tumor volumes measured from 33 CBCT images. The model-calculated and measured tumor volumes, relative to their values on the first RT day, were shown in Fig. 1 With the same optimized parameters, the number of the tumor cells on each of the RT days was calculated using Eqs. (3) and (4), respectively. The computed results were plotted as shown in Fig. 1f. Note that while the Eqs. (3) and (4) were derived with two totally different methods, the numbers of the tumor cells on each day computed by the two equations are almost identical. This has verified that both the analytical and numerical solutions have been implemented correctly in this study.

Impact of tumor doubling time on the derived parameters
In the previous section, the T pot values reported in literature were assigned to the parameter T u . When these values were replaced by the cell division time T div , which is 22 h for adenocarcinoma, the tumor volumes calculated by the kinetic model would be in a zigzag shape as shown in Fig. 2a (in red). The rapid increase after each weekend is due to the short doubling time used in the model, no matter how the other parameters were optimized. This example demonstrated that without using an appropriate tumor volume doubling time, the kinetic model cannot fit to the measured tumor volumes accurately, and T div is not suitable to be used as T u in this model. In contrast, the tumor volume doubling time (T vol = 221.6 days) is much longer than the potential doubling time [39]. Use of the volume doubling time makes the tumor growth over weekend negligible, and the optimized curves are very smooth. However, the calculated volumes still could not fit to the measured tumor volumes (Fig. 2b). Compared to T div and T vol , the measured T pot values showed the best fit for the kinetic model (Fig. 1).
The fitness of different doubling times used in the model was quantitatively assessed for the five patients. It has been found that the T pot -based residual values R are less than or comparable to those generated by the T div or T vol -based model, for each of these patients. All the optimized residual values were shown in Table 3.

Stability of the kinetic model at the optimized parameters
With the proportion of the CSCs (p) changed from 10 − 3 to 1.0, and the other parameters fixed at the optimized values as shown in Table 2, the resultant function R was plotted in Fig. 3a. It can be found that only one minimum point exists for each patient. Similarly, let the probability of the symmetry division (δ) change within a reasonable range of 0.001 to 0.2 as reported in [22,23], and the surviving fraction of CSCs change within the range of 0.001 and 1.0, the objective function R was computed with results illustrated in Fig. 3b and c, respectively. Figure 3 showed that the kinetic model is stable at the optimized points for the parameters p and SF u , where the minimum points for the function R were clearly demonstrated. In contrast, the rate of the CSC self-renewal was shown to be sensitive to the objective function R.

Verification of the cell surviving fraction
Tumor cell radiosensitivity is a major determinant of tumor response to radiation treatment [2]. The following experiment was conducted to verify the surviving fraction of cancer cells under radiation. The cell lines of NCI-H2170, originated from a lung patient with squamous cell carcinoma, were cultured and investigated. The first step of the experiment was to evaluate the growth of this type of cells without irradiation. The cell  Table 4). The second step is similar to step 1 except that the cells were irradiated by 2 Gy after attached to the plate. The numbers of the cells were counted 3 times (days 2, 6, and 9), 3 samples per day.
The ratios of the cell numbers in step 2, relative to those in step 1 are 0.78, 0.56 and 0.44 for the three days.
Due to the effect of cell repopulation, the ratio N 2 /N 1 continuously decreases. Exponential extrapolations showed that the number of the survived cells on day 0 is 1.476 (million) for the non-irradiated experiments, and 1.344 (million) for the irradiated experiments. Consequently, the surviving fraction after 2 Gy, SF 2 , counted as the ratio N 2 / N 1 on day 0, is 0.91. The fitted exponential formula in Fig. 4a showed that the cell division time is ln2/0.2865 = Fig. 1 (a-e) The tumor volumes calculated by the kinetic model compared with those measured from CBCT images for the five patients; (f) The tumor volumes calculated with the analytical and numerical solutions (3) and (4), respectively. Note that the x-axis represents the treatment dates, and the y-axis highlights the range of tumor regression 2.42 days or 58.1 h. The measured values (0.91 and 58.1 h) are consistent to 0.90 for SF 2 and 60 h for the cell division time reported for squamous cell carcinoma [35,36].

Discussion
Radiation dose for locally advanced lung cancer patients is often limited by radiation toxicity to healthy organs due to the large size of tumor volume, and the 5-year surviving rate for these patients is only about 5~14%.
For non-small cell lung cancer (NSCLC) patients, tumor volume usually will shrink after a few treatment fractions. Consequently, the original treatment plan could be reassessed to determine if the remaining tumor can be better targeted [42,43]. The goal, of course, is to use the minimal dose that is necessary to eliminating cancer cells in the remaining tumor target. Cancer stem-like cells complicate treatment planning because of their reported radiation resistance [44]. Kinetic modeling of CSCs may provide a better understanding of tumor response to radiation therapy, and is valuable for developing new treatment regimens for cancer treatment [26].
In the context of CSCs, the minimal dose required may depend on the proportion of CSCs and their radio-sensitivities. The fraction of the stem cells is central to the cancer stem cell debate [26]. Although often a small population has been reported to be present, e.g. 0.2~0.8% for pancreatic cancer [45] and~1% for prostate cancer [46], a much larger fraction (27%) with tumor initiation capability may be present for other cancers, such as melanoma [47]. Indeed, in some tumors almost the entire tumor could be comprised of CSCs. Other factors also could affect the fraction of CSCs in a tumor. For example, CSC numbers could depend on a tumor's stage [19] and the frequency of CSCs might increase during tumor progression [27]. In the current study, the proportion of CSCs derived from the analyses was as high as 40% for the locally advanced NSCLC patients. It was found that the average proportions of stem-like cells were 22.9 and 40.8% for the patients with squamous cell cancer and with adenocarcinoma, respectively. Furthermore, in previous studies it was presumed that the CSC self-renewal rate δ is 0.01 [19,23]. In this study we found that the rate is between 0.018~0.048 for 4 out of the five patients except patient A2 who has a relatively large δ value (0.16). This patient missed treatment for one week after RT started, but the reason for the relatively large δ needs to be further investigated.
While many kinetic models of cancer stem cells have been developed to investigate cell biological dynamics [23,24], quantifying the population of the stem cells, and measuring parameters used in these models have proved to be difficult [21,48]. We extended the kinetic models to include a compartment of dead cells so that changes in tumor volume observed from a sequence of CBCT images can be incorporated into the investigation of radio-biological parameters. This requires the mathematical model to be relatively simple, and in this study, components such as transitional progenitors were not explicitly modeled. For the same reason, many other factors such as cell quiescence, intercellular signaling, immune function and tumor microenvironment were not specifically modeled in this study [49][50][51].  Nevertheless, the set of ordinary differential equations modeled three representative compartments, and characterized the time evolution of each cell subpopulation with the analytic solution aligned with the day-to-day tumor volume changes. As a preliminary study, this work demonstrated the potential to correlate the properties of CSCs with tumor volume changes measured from anatomical images acquired during the course of radiation treatment. It should be mentioned that due to the quality of CBCT images, uncertainties may exist in the contoured tumor volumes, especially for patients with tumor attached to chest walls or mediastinal structures. However, these uncertainties averaged over 33 fractions may have limited impact on the residual function R. Furthermore, with more anatomical and functional images such as magnetic resonance imaging (MRI) and PET/CT being used during the course of radiation therapy, the tumor volumes could be measured more precisely and consequently will help improve the accuracy of the kinetic model. On the other hand, new therapies that target CSCs are under development. Antibody-based constructs that target cell surface proteins expressed on CSCs, natural products like resveratrol and curcumin, some small molecule inhibitors, and classical drugs such as metformin have all shown potential cytotoxicity against CSCs [52]. Since radiation therapy is a standard treatment for many lung cancers, these agents will be adjuvant to RT and since CBCT images are often acquired during the course of fractionated RT, tumor volume data used for the current analysis could be available for subsequent evaluation of the efficacy of new therapies that target CSCs. To our knowledge, this is the first kinetic model that relates the changes in proportion of

Conclusion
A mathematical model has been developed to relate the kinetics of cancer stem-like cell proportions to tumor volume changes observed from a sequence of CBCT images acquired during the course of radiation treatment. The model integrated with clinical data may help in our understanding of the underlying mechanism of tumor response to radiation treatment and perhaps other therapies that target CSCs, and therefore is valuable in the development of new treatment regimens for individual patients.