Normally,cells exposed to ionizing radiation (IR) are blocked at a certain cell cycle phase,which is thought to be caused by the generation of DNA double-strand breaks (DSBs) and forced cells to initiate DNA repair prior to the next round of cell division^{[1, 2]}. G₁ and G₂ arrest following exposure to radiation are widely reported in recent years^{[3, 4]}. In our previous studies,prolonged G₂ arrest was observed besides classical G₁ arrest after X-ray or heavy-ion irradiation of a melanoma cell line 92-1^{[5, 6]}. The underlying mechanism for long-term G₂ arrest is assumed to be the degradation of cyclin B1,an important cell-cycle checkpoint factor. After extended periods of cell cycle arrest,the vast majority of cells eventually become senescent^{[7, 8, 9]}. These findings provide a solid experimental correlation between cell cycle block with cellular radio-sensitivity,which enhances the importance of cell cycle process.

Many models have been developed to describe the dynamical behavior of cell populations^{[10, 11, 12, 13, 14]}. Some of them have the characteristics of good expansibility^[15]. Liu et al^[16] and Dawson et al^[17] constructed a mathematical model to investigate how the stimulus influenced the cell cycle distribution. Brown et al^[18] developed a simulation method to analyze the cell cycle distribution using flow cytometric data. In addition,several population models were used to study the dynamics of the irradiated tumor cells and the various cell fates have been considered^{[19, 20, 21, 22]}. However,there are few studies focusing on the cell cycle progress following ionizing radiation and the relevant experimental data are lacking. In the present study,a computational model was developed by designing a phenomenological system,which is consistent with our experimental data. The results should provide a first step to more precise prediction of the cell status following exposure to ionizing radiation^{[23, 24]}.

Human choroidal malignant melanoma 92-1 cells were cultured in RPMI 1640 medium (Sigma,USA) complemented with 10% fetal bovine serum (FBS,USA),100 g/ml streptomycin,and 100 U/ml penicillin in a humidified atmosphere of 5% CO₂. Cells were seeded in 12-well plates 2 d before irradiation,and grown to <70% confluence at the time of irradiation. Cells were exposed to X-rays (100 kV,5 mA,no additional filter; Faxitron RX650,USA) at room temperature and a dose rate of 1 Gy/min.

For cell synchronization,92-1 cells were grown in medium containing 1 mmol/L thymidine (Sigma,USA) for 24 h without delay irradiation.

Cells were harvested and fixed with 70% pre-chilled ethanol for at least 24 h at -20 ℃. Prior to analysis,fixed cells were washed twice with PBS,centrifuged at 300×g for 5 min,treated with 100 g/ml. RNase A and 50 g/ml propidium iodide (BD Biosciences,UAS) for 30 min at 37℃,and analyzed by FACS calibur flow cytometry (Becton Dickinson,USA). DNA content was measured with CellQuest (Becton Dickinson,USA) and cell cycle distribution was analyzed by FlowJo software.

A cell-cycle model with three phases is proposed,G₀/G₁ phase,S phase,G₂ phase and M phase. The duration of each phase is assumed to obey Gaussian distribution. In our model,there were N variables S_n (n=1…N) represented the cell cycle phase of N different cells. Their value S_n=0,1,2,3means that the cell n was in S,G₂,M or G₁ phase,respectively. For the n^th cell,four random numbers (t_n⁰,t_n¹,t_n²,t_n³) are created to describe the duration of S,G₂,M and G₁ phase of the cell. These random numbers obey fixed Gaussian distribution. Because of the efficiency of synchronization in the actual experiments was not ideal,the synchronized cells were not actually synchronized at exactly the same time-point in cell cycle duration even though their cell cycle measured with flow cytometry seemed to be at the same phase. Thus,another random number t_n is created to represent the time-point in the cell cycle of cell n in simulations. According to the experimental data (the proportions of cell phases at time point 0 h),the time-point t_n of one cell in the cell cycle prior to irradiation is also supposed to obey Gaussian distribution.

DSBs induced by ionizing radiation have to be widely repaired^[25]. In the simulation,the numbers of DSBs and other DNA lesions for each cell should obey a Poisson distribution and the mean of this distribution should be proportional to the radiation intensity^[26]. If the repair of same DNA lesions costs cell the same time to complete,the total recovery time for one cell t^repair will obey a Poisson distribution,too. It is demonstrable that lager dose radiation could cost cells more time to repair. Although the kinetics of DSB repair follows a double exponential curve^[27] and seems to be independent of radiation intensity after exposed to low-dose radiation^[28],the repair time of the cells pretreated with drug after exposure to low-dose radiation seems to be proportional to radiation intensity reversely^[29]. Moreover,besides the DSB,the kinetics of other DNA lesions repair follows a linear curve^[30]. The relation between repair time and radiation intensity is still not clear^[31],especially after the moderate-dose and high-dose radiations. Here we assumed that expectation of t^repair should be proportional to radiation intensity (t^repair∝dose) and the average recovery time was supposed as t^repair(h)=2×radiation intensity(dose). Thus,when cell n suffers radiation,a random number t_n^repair,which obeys a poisson distribution,will be created to describe the duration of DNA repair. Based on our previous results,melanoma 92-1 cells are blocked in G₂ phase by ionizing radiation^[9]. The mean and variance of the distribution are proportional to the irradiation dose. For one cell,two unique thresholds (Tr₁,Tr₂) are defined: if t_n^repair<Tr₁,the cell will completely repair the radiation-induced damage; if Tr₁<t_n^repair<Tr₂,the cell can pass the G₂ checkpoint but will be blocked at the following G₁ checkpoint; and if Tr₂<t_n^repair,there will be a prolonged G₂ arrest for this cell. The G₁ checkpoint is set at 2 h before the end of G₁ phase. A completely repaired cell will come back to normal cell cycle and have a division capacity at M phase. In this model,when the cell divides into two cells,the daughter cells will have the same duration of each phase as mother cell.

Insimulations N=300,and parameters are shown in Table 1. For cell n,t_n in simulation correspond to t_n h. One time step was set as 0.01 h. That is,in the simulation process of our model we let t_n+0.01→t_n after one time step to present growth of cell n. Totally 10 000 time steps are taken for simulation. Figure 1 showed the pseudo code and detailed simulation rules.

表 1(Table 1)

Table 1 Parameters used in our simulations. Average durations of cell-cycle phases were empirical and the variations of them were assumed. Two thresholds were assumed based on the experimental data

Table 1 Parameters used in our simulations. Average durations of cell-cycle phases were empirical and the variations of them were assumed. Two thresholds were assumed based on the experimental data

图1

Fig. 1

Figure 1 Pseudo code of model

In order to give a qualitative explanation for the cell cycle dynamics,a theoretical model is constructed to explain the long-time asymptotic behaviors of cells that suffered from IR. If the cell death is neglected,three ordinary differential equations (ODEs) are written to describe the dynamics of cell cycle as following^[22]:

In order to obtain the cell cycle dynamics of the cells in difference phases,92-1 cells were synchronized in late G₁ phase and exposed to various doses of X-rays immediately after release and then their cell cycle distribution was monitored with flow cytometry up to 96 h.

Although 80% cells were synchronized at the later G₁ phase,only 10% cells remained in G₁ phase 6 h later after release no matter that cells were exposed to ionizing radiation or not (Figure 2). Within the time range of the experiments (0-96 h),the proportion G₀/G₁ phase cells in the whole cell populations at normal condition displayed an oscillation at 48 h after release from thymidine treatment and then reached to a stable phase. 24 h later after exposure to 2 Gy of radiation,about 80% cells returned to G₁ phase and remained unchanged then. Upon exposure to 5 Gy of radiation,only about 30% cells returned to G₁ phase. No cells returned to G₁ phase after exposed to 10 and 15 Gy of radiation.

图 2(Figure 2)

Fig. 2

Figure 2 Dynamics of cell cycle distributions at different doses A. Dynamics of G1 phase distributions; B. Dynamics of S phase distributions; C. Dynamics of G2/M phase distributions

No matter irradiated or not,all cell populations exhibited a peak in S phase 6 h after cell cycle release (Figure 2). Cells without irradiation displayed a second peak at 36 h after release,while those cells who were exposed to radiation displayed a low proportion of S phase cells without any obvious peak.

As shown in Figure 2,no cells were found in G₂/M phase after synchronization. However,the proportion of G₂/M phase cells increased right after release. For control cells,G₂/M phase cells reached to a peak of 35% 6 h after release and then decreased to less than 10%. However,it increased again 24 h later and reached to another peak at 36 h after release.After 2 Gy irradiation cells,the proportion of G₂/M phase cells reached to a peak of 50% and then decreased gradually to less than 10%. For 5 Gy irradiation,the proportion of G₂/M phase cells reached to a higher peak but decreased slightly. After 10 or 15 Gy irradiation cells,the proportion of G₂/M phase cells reached to a peak as high as 75% and remained no changed during subsequent culture.

Obviously,cells without irradiation entered S phase and G₂/M phase right after release. X-ray radiation blocked cells in G₂/M phase in a dose-dependent manner. If the irradiation dose was low,a small percentage of cells returned to the G₁ phase. However,when the dose was over 10 Gy,the majority of the cells were arrested in a prolonged G₂ phase. These results were consistent with our previous studies.

In the simulations,the comparisons between simulation results and experimental data were shown in Figure 3. The degree of fitting was calculated to test the effectiveness of this model. Table 2 showed that there was a consistency between simulation results and experimental data at 0 and 2 Gy. As the irradiation doses were increased,the proportion of S phase cells between simulation results and experimental data showed difference and the amount of proportion of S phase cells was about 10% after 24 h while in simulation results,it closed to zero. It implied that a slight S phase arrest might exist in high-dose irradiated cells,which had not been considered in our model.

图 3(Figure 3)

Fig. 3

Figure 3 Comparison between simulation results and experimental data at different doses A. 0 Gy; B. 2 Gy; C. 5 Gy; D. 10 Gy

表 2(Table 2)

Table 2 The degree of fitting between simulation results and experimental data for different cell phases at different doses,the formula used here is R2=1－|SST－SSE|/SST,SST was total sum of squares,SSE was error sum of squares

Table 2 The degree of fitting between simulation results and experimental data for different cell phases at different doses,the formula used here is R2=1－|SST－SSE|/SST,SST was total sum of squares,SSE was error sum of squares

The influence of parameters (Tr₁,Tr₂) was studied on the proportion of prolonged G₂-arrested cells after radiation. The variations of the proportion of prolonged G₂-arrested cells with different Tr₁ and radiation doses were given when the cells are supposed to be ideally synchronized between G₁ and S phases (Figure 4). Tr₁ was connected with cellular radio-sensitivity of cells. If cells had a higher Tr₁,the proportion of prolonged G₂-arrested cells would be 0. If the radiation dose was increased,there would be more prolonged arrested cells and the proportion would reach to 1. 92-1 cells are very radiosensitive,so Tr₁ for them is small. In our model the value was given as 1 h.

图 4(Figure 4)

Fig. 4

Figure 4 Influences of parameters in phenomenological model and theoretical model A. Heat map of proportion of prolonged G2-arrested cells with different Tr1 and radiation intensity; B. Heat map of ln(G01/G2M) with different parameters a and c

From the Eq.(1-3),the asymptotic behavior of G₀₁/G_2M when time t tends to infinity (t→∞) can be analytically given by

Then the values of G₀₁/G_2M was estimated at 96 h after radiation based on the experimental data. From the heat map,the estimated values of our parameters a and c at different doses of radiation are obtained.

[LL] When dose was null,the parameters a and c were obtained from normal cell cycle. When dose was 2 Gy,parameter c did not change while a dropped to a small number (see in Table 3),which means a strong G₁ arrest. When the dose ascended,parameter c was falling,which means a strong G₂ arrest.

表 3(Table 3)

Table 3 Parameters in the theoretical model

Table 3 Parameters in the theoretical model

Several mathematical tools were designed to depict the dynamics of cell phase curves^{[32, 33, 34]}. Enderling et al^[35] had developed a model for tumor survival after radiotherapy. Ochab-Marcinek et al^[36] reported a model to discuss the influence of radiations on cell-cycle delay. But a theoretical model for cells containing cell-cycle arrest and DNA repair process has not been established yet. The dynamics of radiation-induced cell-cycle arrest and DNA repair process was obscure and hard to use a group of simple differential equations to be described. So,it was difficult to get analytic solutions for these cell behaviors^[37]. A phenomenological simulation is a good alternative tool. As the results shown in this study,the simulations semi-quantitatively fit experimental data if proper parameters were chosen. Because of tumor cell lines vary in radiosensitivity and DNA repair capacity,our next goal is to verify its feasibility of this model in describing the dynamics of other types of human cell lines.

In our previous studies,different melanoma cell lines showed different radio-tolerance and response after IR-irradiation . In this study,92-1 cell line was chosen,which showed a prolonged G₂ block after a high dose ionizing radiation. The proportions of every cell phase at several time-points were measured. Some cells pass G₂ arrest,but they were still blocked at G₁ checkpoint,especially when the dose was less than 5 Gy. G₁ checkpoint was more stringent than G₂ checkpoint^[38]. Cells undergoing incomplete DNA damage repair might pass the G₂ checkpoint and be blocked again at G₁ checkpoint,thus express a cell cycle suspension. This phenomenon was considered in our models. Our simulation results were consistent with the experimental data would support this point in reverse. In other papers,cells treated with drugs could be stopped at G₁ or G₂ checkpoints with a probability closed to 1^[23],which is similar in proportion to high-intensity irradiation in our experiments.

Although synchronized cells were used in this work,this method could also be extended to simulate cell populations in any distribution of cell cycle. If there were other repair modes such as mismatch repair in cells,the duration of repair must be rewritten to a complex mode. Several publications had shown that there are differences between ionizing radiation and ultraviolet radiation from DNA repair mode to molecular response^{[39, 40]}. Next step,we will develop our phenomenological model to describe cell behaviors after exposure to other radiation types such as ultraviolet exposure.

A phenomenological model was constructed in this paper,which could be used to describe the dynamics of cell cycle of cells exposed to ionizing radiation and was supported by the experimental data. Our model could be used to predict dynamics of cell cycle of different cells exposed to a certain irradiation dose. Because this model is easy to run by written code,it has a good expansibility for studying the behaviors of cell populations under various conditions.