Cell-Cycle Synchronization Prior to Radiotherapy: A Mathematical Model of the Use of Gemcitabine on Melanoma Xenografts

Abstract

Radiotherapy can differentially affect the phases of the cell cycle, possibly enhancing suppression of tumor growth, if cells are synchronized in a specific phase. A model is designed to replicate experiments that synchronize cells in the S phase using gemcitabine before radiation at various doses, with the goal of quantifying this effect. The model is used to simulate a clinical trial with a cohort of 100 individuals receiving only radiation and another cohort of 100 individuals receiving radiation after cell synchronization. The simulations offered in this study support the statement that, at suitably high levels of radiation, synchronizing melanoma cells with gemcitabine before treatment substantially reduces the final tumor size. The improvement is statistically significant, and the effect size is noticeable, with the near suppression of growth at 8 Gray and 92% synchronization.

1. Introduction

Melanoma is a common form of skin cancer and has the highest mortality rate of skin cancers in the U.S. [1]. Radiation therapy works by inducing cellular damage and thus has the potential to damage normal tissue [2]. Historically, melanoma has been deemed a radioresistant tumor [3,4]. Radiation is used when surgery is not possible and sometimes in combination with other therapies [5,6].

To proliferate, cancer cells progress through the cell cycle, some stages of which are known to be more sensitive to radiotherapy [7,8,9,10,11,12,13]. A variety of drugs are known to produce cell-cycle arrest in tumors, resulting in the synchronization of the cells in a particular stage of the cell cycle [14,15,16,17]. By administering a synchronizing drug to increase the proportion of tumor cells in a radiosensitive stage of the cell cycle, more cells may die from radiotherapy than in an unsynchronized tumor. As such, cell-cycle synchronization prior to radiation treatment is regarded as a potential therapeutic target, which has been explored in several studies [18,19].

For example, cell survival from radiotherapy was lowest during mitosis and late $G_1$ to early S phase in HeLa cells, a cervical cancer cell line [20]. Otani et al. carried out a similar study, exposing various human cell lines, including HeLa (cervical cancer), HCT166 (colorectal cancer), PANC-1 and MIA PaCa-2 (pancreatic cancers), MG-63 (osteosarcoma), and B16BL6 (melanoma), to fixed doses of radiation [19]. PANC-1 and MG-63 both showed some increased susceptibility to radiotherapy during the S phase, but this was not statistically significant. B16BL6 was uniquely radioresistant among the cell lines tested and at high doses of radiotherapy showed significant susceptibility in the S phase relative to other cell-cycle phases. The remaining cell lines were highly radiosensitive, with all cells dying at moderate doses of radiation, and they did not show clear susceptibility in a given cell cycle stage.

An observed improvement in radiation response, even if statistically significant, may not make a big enough difference to be worth pursuing in a clinical trial. Predicting a likely outcome requires a quantitative approach. Some attempts have been made to quantify aspects of this phenomenon using mathematical models [21,22,23].

The mathematical model used here to investigate cell-cycle effects on radiosensitivity is constructed in a nested manner to take advantage of monolayer, spheroid, and xenograft data, based on the work of He et al. [24]. In some studies, these three experimental modalities are investigated simultaneously, although this is rare [25]. Instead, the mathematical model here relies on a series of experiments on melanoma cell lines reported in the literature.

A three-dimensional, or spheroid, cell culture is an important intermediate step between monolayer and xenograft models [26,27,28,29,30,31,32]. To represent a spheroid culture, the model used in this study retains the parameters for the cell cycle determined from the monolayer experimental data and adds quiescent and necrotic compartments. In addition, the model incorporates the cytokine, tumor necrosis factor alpha (TNF-$\alpha$), which is produced by the necrotic core and induces apoptosis in proliferating cells as described by Wallace et al. and He et al. [24,33]. Spheroid growth data are taken from Browning et al. [34].

The xenograft model retains the features and parameters of the spheroid model, adding in angiogenesis and its dual roles of providing nutrients and removing cellular waste. Vasculature grows in the model in response to vascular endothelial growth factor (VEGF), produced by quiescent and proliferating cells in response to TNF-$\alpha$. In addition to delivering nutrients, vasculature removes waste and clears both TNF-$\alpha$ and VEGF. The full model is explained in detail by He et al. [24]. Xenograft control data are taken from Dorasamy et al. [35].

Radiation treatment is modeled after the experiments of Otani et al., in which gemcitabine, a clinically used anti-cancer drug, induces cell-cycle arrest during the early S phase in the B16BL6 melanoma cell line [19]. After cell-cycle arrest, various doses of radiotherapy were administered, as measured in Grays (Gy), the standard unit of ionizing radiation. An in vitro analysis showed higher radiosensitivity of cells in early S-phase as compared to the other cell cycle stages at a dose of 8 Gy. The in vivo study compared a single dose of radiation given to B16BL6 xenografts in mice with and without prior gemcitabine treatment. This experiment resulted in regrowth of a smaller tumor in the gemcitabine-synchronized xenograft tumors.

One clinical study compared the radiation treatments of 32 Gy in 8 Gy weekly fractions versus 50 Gy in 2.5 Gy daily fractions for melanoma. Not much difference in response between these dosage patterns was observed [36]. For this reason, and because 8 Gy produced more noticeable results in vitro, our numerical experiment administers 24 Gy in 8 Gy weekly fractions for a total of three treatments applied to the xenograft model, with and without the reported synchronization from gemcitabine.

2. Methods

Unlike most modeling approaches that seek to parameterize a model using a single data set, a nested set of models is matched to monolayer, spheroid, and xenograft data in that order. In doing so, the model leverages more experimental data and reduces the number of parameters that need to be fit at each level of the model. The linear monolayer model represents proliferating cells transitioning through the cell cycle, and it grows exponentially. This system is expanded to a non-linear spheroid model by including quiescent and necrotic cells, as well as the production and action of TNF-$\alpha$ via Michaelis–Menten kinetics, giving bounded growth. The spheroid model is further expanded to the xenograft model with the inclusion of growing vasculature through angiogenesis. It also incorporates the production of vascular endothelial growth factor (VEGF) by proliferating tumor cells interacting with TNF-$\alpha$ via the p75 receptor and by quiescent cells presumed to be hypoxic. The full xenograft model grows quickly at first, then settles into a linear growth pattern. The equations for the system are similar to those of He et al. [24] and may be found in Appendix A; a more thorough description of the parameter selection and the development of the model equations can also be found in the paper by He et al. A compartment diagram of the system is shown in Figure 1.

Monolayer data, including doubling time and the proportion of cells in each phase of the cell cycle, are used to determine the transition parameters for the proliferating cells. These parameters are retained in a spheroid model, which is further matched to determine transition parameters for quiescent and necrotic cells, as well as TNF-$\alpha$ production and interaction with proliferating cells. The spheroid parameters are finally lifted to a xenograft model, which is matched with tumor growth data to give the remaining parameters. Appendix B shows the parameters used for the untreated xenograft model and both submodels (

Table A1. Model Parameters.

	Monolayer	Spheroid	Xenograft	Explanation
$d_a$	0	0	0	Natural Death Rate
$c_1$	1.53372	1.53372	1.53372	$G_1$ to S Transition
$c_s$	1.5004	1.5004	1.5004	S to $G_2$ Transition
$c_2$	63.104	63.104	63.104	$G_2$ to $G_1$ Transition
$v_0$	N/A	2.9	100	Ambient Nutrient
j	N/A	1.1	1.1	TNF-$\alpha$ Production
k	N/A	0.14	0.14	TNF-$\alpha$ Removal
$s_1$	N/A	1	1	B Half-Saturation Constant
$s_q$	N/A	300	300	C Half-Saturation Constant
e	N/A	0.08	0.08	Necrotic Death of Q
m	N/A	0.004	0.004	Natural Removal of N
$s_n$	N/A	60	60	F Half-Saturation Constant
$c_q$	N/A	55	55	Maximum Rate of C
$c_f$	N/A	0.83	0.83	TNF-$\alpha$-Induced Apoptosis of P
$c_v$	N/A	N/A	70	Vascular Growth
$c_R$	N/A	N/A	223.05148	Production of VEGF by P
$s_R$	N/A	N/A	103.45178	Stimulation of VEGF by TNF-$\alpha$
$q_R$	N/A	N/A	5.50064	VEGF Removal
$c_{qv}$	N/A	N/A	0.01	Production of VEGF by Q
$c_{qs}$	N/A	N/A	0.17877	Controls C
$s_v$	N/A	N/A	100	Controls Vasculature
$s_h$	N/A	N/A	0.1	H Half-Saturation Constant
$f_h$	N/A	N/A	0.1	Controls H
$g_h$	N/A	N/A	0.1	Controls H

), as well as the initial conditions (

Table A2. Initial Conditions.

	Monolayer	Spheroid	Xenograft	Explanation
${G1}_0$	$2.8735\times {10}^5$ cells	0.00908 mm3	69.76 mm3	$G_1$ Phase Cells
$S_0$	$2.08\times {10}^5$ cells	0.00657 mm3	50.50 mm3	S Phase Cells
${G2}_0$	$0.049\times {10}^5$ cells	0.000155 mm3	1.190 mm3	$G_2$ Phase Cells
$Q_0$	N/A	0	0	Quiescent Cells (mm3)
$N_0$	N/A	0	0	Necrosis (mm3)
$T_0$	N/A	0	0	TNF-$\alpha$
$R_0$	N/A	N/A	0	VEGF
$V_0$	N/A	N/A	0	Vasculature

).

2.1. Monolayer Parameters

Linear models of the monolayer culture show spontaneous desynchronization of the cell cycle [33,37]. This is because the limiting cell-cycle distribution is the eigenvector associated to the largest eigenvalue of the system, which is given by the intrinsic growth rate, computed from the doubling time as derived in Wallace et al. [33]. This approach allows empirically derived measures to be used to calculate all transition rates for the linear monolayer model. The measurements were determined from three studies—the doubling time was taken from Liberato et al. [38], the cell-cycle analysis from Pan et al. [14], and the initial conditions from Otani et al. [19].

2.2. Spheroid Parameters

Spheroid parameters were derived from a fit of data from Browning et al. [34]. They took confocal images of 318 WM938b spheroids and 317 WM793b spheroids, both of which are human melanoma cell lines. The WM793b spheroids did not reach a steady state during their experimental window but the WM938b spheroids did, so the second set was selected. Additionally, they studied three different seeding densities for their spheroids (2500 cells, 5000 cells, and 10,000 cells); the dense spheroids were close to a steady state in the plots they presented, while the 2500 cell spheroids started with a smaller volume and moved toward a steady state throughout the experimental window and were thus the better choice for parameterizing the spheroid model.

The experimental measurements include R (total spheroid radius), $\varphi$ (fraction of the radius from the center of the spheroid to the outer edge of the quiescent compartment), and $\eta$ (fraction of the radius accounted for by the necrotic core). The authors provided a link to their github that includes the raw measurements of the spheroids from confocal images on experimental days 3, 4, 5, 7, 10, 12, 14, 16, 18, 21 [34]. Volumes for the whole spheroid (W), proliferative compartment (P), quiescent compartment (Q), and necrotic compartment (N) were computed, assuming that the spheroid is a perfect sphere comprised of a necrotic core, covered by quiescent cells, which are in turn covered by proliferating cells. The resulting volumes were divided by ${10}^9$ so they would be in units of mm³ and are reported in Appendix C (

Table A3. Experimental spheroid volumes (mm³) [34].

Day	${\mathit{W}}_{\mathit{vol}}$	${\mathit{P}}_{\mathit{vol}}$	${\mathit{Q}}_{\mathit{vol}}$	${\mathit{N}}_{\mathit{vol}}$
3	0.0158	0.0158	0.0001	0
4	0.0215	0.0212	0.0002	0
5	0.0262	0.0258	0.0004	0
7	0.0477	0.0406	0.0071	0
10	0.0903	0.0428	0.0434	0.0041
12	0.0922	0.0395	0.0436	0.0091
14	0.1287	0.0479	0.0617	0.0192
16	0.1475	0.0528	0.0638	0.0308
18	0.1503	0.0471	0.0605	0.0427
21	0.1669	0.0448	0.0540	0.0680

). Mean volumes across the spheroids at each time point were computed, and the model parameters were manually tuned to match the data for whole volume, proliferating, quiescent, and necrotic volumes simultaneously. The match to data is shown in Figure 2A.

2.3. Xenograft Parameters

The remaining parameters for the xenograft model without radiotherapy were determined from a fit to the data in Dorasamy et al. [35]. They used A375 xenografts—a different human melanoma cell line. They stated that after 14 days of tumor growth, the xenografts attained a volume of $1090.0\pm 134.7$ mm³; the MATLAB “grabit ” function was used to extract data for the remaining experimental timepoints [39]. The extracted values for mean xenograft volume and standard error of the mean are presented in Appendix C (

Table A4. Experimental in vivo volumes (mm³) [35].

Day	Volume (Mean ± SEM)
1	$121.4$
4	$227.7\pm 58.8$
6	$397.3\pm 52.7$
8	$538.9\pm 83.1$
11	$771.6\pm 99.8$
14	$1090.8\pm 134.8$

).

It is of note that all but one of the spheroid parameters was held constant in parametrizing the xenograft; $v_0$ was allowed to vary because the nutrient availability in the microenvironment is expected to differ when transitioning from the cell culture to an in vivo model. Parameters were first fit manually and later refined using a genetic algorithm with 200 individuals over 200 generations. The relative mean squared error (rMSE) of the manual parameterization is 0.076709; the rMSE after the genetic algorithm is 0.03773. The plot of the final data fit for the xenograft is shown in Figure 2B. After the genetic algorithm, the model curve for the whole tumor volume fits within the SEM boundaries for each time point.

2.4. Numerical Treatment Simulations

There are two elements to the treatments administered by Otani et al.: cell-cycle synchronization in the S-phase by the administration of gemcitabine and irradiation with varied Grays (Gy) of ionizing radiation. Radiotherapy is simulated at 2, 4, and 8 Gy according to the in vitro cell survival measurements for the various compartments, extracted from the figures in the work of Otani et al. using MATLAB’s grabit function [39]. The cell-cycle phases presented by Otani et al. are not exact—late G0/G1 is taken as $G_1$ in the model, retained G0/G1 as Q, early S as S, and late S/G2/M as $G_2$.

Table A5. Monolayer surviving fraction [19].

Grays	W	${\mathit{G}}_1$	Q	S	${\mathit{G}}_2$
0	1	1	1	1	1
2	0.9896	0.9722	0.8403	0.9810	1
4	0.6643	0.7573	0.6314	0.5037	0.8007
8	0.1964	0.3175	0.2449	0.03154	0.1804

in Appendix C displays the extracted values for each cellular compartment for the different irradiation levels. For the simulations, cells that die due to radiation are immediately added to the necrotic compartment, and surviving cells are reallocated based on the ratios described in Table A5. Per the experimental data from Kaeppler et al., even under 15 Gy of radiation, the vascular structure in the tumors was minimally disrupted, even with the death of endothelial cells [40]. For this reason, the vasculature compartment (V) is not reduced during radiotherapy in the model.

The experimental protocol from Otani et al. involved synchronizing the cells via an intraperitoneal injection of gemcitabine. The cells stayed synchronized for 12–18 hours, after which point they rapidly desynchronized. Mice were separated into two groups: an experimental group that received gemcitabine, followed by irradiation 12–15 hours later (when the cells were maximally synchronized), and a control group that received radiotherapy first, followed by gemcitabine 8–10 hours later to control for the possible effects of the drug. For the in silico model, the experimental group was synchronized and immediately “treated” with radiotherapy, while the control group first received radiation, followed by immediate synchronization of the cells in S-phase.

Figure 3d from the paper by Otani et al. gives us information about the effectiveness of the gemcitabine treatment in synchronizing cells in vivo. Two-photon microscopy was used to detect green fluorescing Fucci-expressing B16BL6 cells, which correspond to the early S-phase when the cells had 2N DNA content. From 12 to 18 hours post-gemcitabine treatment, all viable cells were roughly 92% synchronized in the early S-phase as determined by extracting data corresponding to the mean at 12 hours using MATLAB’s grabit function [39].

Otani et al. did not specify how many Grays the mice received, only stating that radiotherapy was administered at a rate of 0.7 Gy/min, without stating the duration of exposure. Radiation was administered 7, 14, and 21 days after tumor induction. For the purpose of the model, the tumor size at model day 0 corresponds to the size of a tumor one week after induction, so the administration of the treatment on model day 0 is the same as treating on experimental day 7. MATLAB’s grabit function was used to extract the resulting data mean and error bars for the synchronization–irradiation experiment [39].

Simulations of the experiment were run with synchronization at 92% of all viable cells ($G_1$, S, $G_2$, and Q)—92% of the total cell count were moved to the S phase, with the remainder allocated to the other viable compartments based on the ratios of cells in these compartments just prior to synchronization. At 4 Gy, the curves from the model matched the data from the paper satisfactorily. Note that the starting volume was retained from the xenograft fit from Dorasamy et al., as it fell near the datapoint for the tumor volume on day 0 (one week post-tumor induction) presented by Otani et al. [19,35]. The model’s output for 92% synchronization and 4 Gy administered on days 0, 7, and 14 is shown along with the extracted data points in Figure 3.

Different levels of synchronization after treatment with gemcitabine may be present in other contexts, such as human tumors. The model was tested for various levels of synchronization at 0, 2, 4, and 8 Gy, with the final tumor size and percent of control shown in Figure 4.

With the model now parameterized, it is necessary to compare treatment with 92% synchronization at 8 Gy with treatment at 8 Gy with no gemcitabine at all as would be compared in a clinical trial. A simulated clinical trial was conducted by varying parameters in the underlying model to create 200 individuals to be assigned randomly for experimental and control treatment. One hundred of these individuals were treated at 8 Gy with 92% synchronization. The other 100 individuals were treated at 8 Gy with no synchronization (no gemcitabine). Patient individuality was simulated by varying model parameters within a 20% box about the original parameter value, using a Latin hypercube sampling system. Each simulated patient had the same initial conditions (i.e., same initial tumor size), but the variable model parameters resulted in varying the tumor growth over time. The results of this simulation with 100 patients per group are shown in Figure 5.

The resulting final tumor sizes were compared and are illustrated as a histogram in Figure 6.

3. Results

Figure 2A shows the result of matching the underlying spheroid and xenograft models to data from two studies. The spheroid model is matched to measurements by Browning et al. of the proliferating, quiescent, and necrotic compartments, as well as the total spheroid volume [34]. The model is extended to include vasculature, resulting in the xenograft model, shown with data points for the total tumor volume from Dorasamy et al. [35] in Figure 2B.

Figure 3 shows the model comparison for 4 Gy radiation with gemcitabine. The ”unsynchronized” group received gemcitabine after radiation, while the ”synchronized” group received gemcitabine prior to radiation. Data points and error bars are from Otani et al. [19]. Both treatment conditions fail to control the tumors and thus correspond to tumor growth throughout the 21-day simulation, but the synchronized treatment is more effective, resulting in slower tumor growth.

The result of varying both the radiation dose and the level of synchronization is shown as a heat map in Figure 4A. As the level of synchronization increases, the final tumor volume rises for the four control groups. This is likely because the simulation transfers cells from the quiescent compartment into the proliferating S-phase of the cell cycle. This design choice is based on the experimenters’ statement that 92% of all viable cells were in the S-phase following synchronization. The same phenomenon happens for the experimental group at 0, 2, and 4 Gy. At 8 Gy, however, higher synchronization in the experimental group leads to a lower final tumor volume. For 4 and 8 Gy, the second heat map (Figure 4B). shows the ratio of the final tumor volumes, experiment divided by control. The experiment simulation (synchronized) always corresponds to slower tumor growth than the control simulation (unsynchronized), with the best effect at 8 Gy and 100% synchronization.

An in silico clinical trial was carried out with 100 individuals assigned to the control (24 Gy radiation fractionated once weekly in 8 Gy doses) and 100 individuals assigned to the experimental condition (the same radiotherapy following 92% synchronization in the S-phase with gemcitabine). Tumor growth was suppressed in the experimental condition, as shown in Figure 5. The simulated clinical trial shows a clear difference in the outcomes for synchronized and unsynchronized treatments, seen in Figure 6. Synchronized cells had a mean final volume of 104.92 mm³, (sd = 15.45). Unsynchronized cells had a mean final volume of 466.31 mm³, (sd = 58.23). The two means are significantly different ($p<0.0001$).

4. Discussion

As seen in Figure 2, the model presented in this study matches the data from spheroid and xenograft experiments, building on the previous monolayer and spheroid model set forth by Wallace et al. [33]. Similar to the work of Wallace et al., the model is constructed in tiers. In each tier—monolayer, spheroid, and xenograft—additional parameters are incorporated that express the increasing complexity of the tumor environment.

The model is robust across a variety of growth scenarios. Figure 3 shows the fit of the model to experimental data under two simulated treatments: treatment with 4 Gy of radiotherapy preceded by 92% cell synchronization in the S-phase, as well as treatment with 4 Gy of radiotherapy followed by cell synchronization. The model parameters used in Figure 3 are identical to those used in Figure 2B.

The model’s versatility is made possible by its nested construction. Its three functional tiers—monolayer, spheroid, and xenograft—each layer on new model parameters, and each successive tier retains the parameters of the previous tiers. In this way, the complete 24-parameter xenograft model only involves fitting 11 parameters—the remaining 4 that were fit to the monolayer data, and the 9 that were fit to the spheroid data, are fixed. This reduces the potential for over-fitting to a single dataset.

The simulation at 4 Gy and 92% synchronization is a good match to the data in the work of Otani et al. [19]. Figure 3 shows that the 4 Gy dose is not effective at containing the tumor size, either in simulation or in the data. Humans can handle doses of 8 Gy, however, and the simulation makes the advantages of the higher dose clear. Figure 4 simulates an 8 Gy weekly treatment, which is highly effective at reducing the tumor volume when compared to the lower Gy doses. In addition, the simulations show a reduction in tumor volume when cells are synchronized before treatment.

Otani et al. reported a statistically significant difference between single versus combined radiation/synchronization therapies in suppressing the proliferation of cells in vitro. They further reported a statistically significant difference in tumor growth in vivo ($p<0.05$) for their small sample of xenografts. The in silico clinical trial in this study (see Figure 5) confirms the significance of this treatment schedule at 8 Gy ($p<0.0001$) with a larger sample of 200 simulated patients and suppression of tumor growth in the experimental group.

Table A5 in Appendix C shows the highest death rates across all phases of the cell cycle at the 8 Gy dose. Likewise, Figure 5 simulates this 8 Gy dose over 3 weeks for a tumor that is treated without synchronization versus with synchronization. The simulations confirm the results in Otani et al. [19]—a treatment schedule that includes cell synchronization before radiotherapy can help to reduce tumor growth.

Future Directions

Otani et al. found that up to 92% of all xenograft cells were found in the S-phase after gemcitabine treatment; this was only possible in our model if we allowed quiescent cells that reside in the $G_0$ phase described by Otani to also enter the S-phase upon synchronization. This assumption should be validated experimentally with further cell-sorting assays to determine whether quiescent cells reenter the cell cycle and arrest in the S-phase or if they are simply not responsive to gemcitabine. The model could easily be adjusted to reflect a different assumption if warranted by experimental data. The experimenters also claim that gemcitabine alone suppresses growth in vitro, although the exact method of measuring this was not given.

Additionally, a xenograft differs substantially from a human tumor, as it lacks an immune response and is induced in a mouse, which has a different metabolism and size from a human. Future work would add a simplified immune response to the system of differential equations underlying this simulation. At present, how to measure immune system parameters directly remains a quandary. The immune response would naturally enhance the effect of radiation and could possibly show tumor suppression at a lower dose of radiation, or eradication at a higher dose.

5. Conclusions

Is it worth synchronizing cells before administering radiation treatment? The simulations offered in this study support the statement that, at suitably high doses of radiation, synchronizing melanoma cells with gemcitabine before treatment results in substantial improvement in the final tumor size. The improvement is statistically significant and the effect size is noticeable, with near suppression of growth at 8 Gy and 92% synchronization.