- Research Article
- Open Access
Cellular automaton-based model for radiation-induced bystander effects
BMC Systems Biology volume 9, Article number: 90 (2015)
The radiation-induced bystander effect is a biological response observed in non-irradiated cells surrounding an irradiated cell. The bystander effect is known to be induced by two intercellular signaling pathways, the medium-mediated pathway (MDP) and the gap junctional pathway (GJP). To investigate the relative contribution of each signaling pathway, we have developed a mathematical model of the cellular response through these two pathways, with a particular focus on cell-cycle modification.
The model is based on a cellular automaton and consists of four components: (1) irradiation, (2) generation and diffusion of intercellular signals, (3) induction of DNA double-strand breaks (DSBs), and (4) cell-cycle modification or cell death. The intercellular signals are generated in and released from irradiated cells. The signals through the MDP and the GJP are modeled independently based on diffusion equations. The irradiation and both signals raise the number of DSBs, which determines transitions of cellular states, such as cell-cycle arrest or cell death.
Our model reproduced fairly well previously reported experimental data on the number of DSBs and cell survival curves. We examined how radiation dose and intercellular signaling dynamically affect the cell cycle. The analysis of model dynamics for the bystander cells revealed that the number of arrested cells did not increase linearly with dose. Arrested cells were more efficiently accumulated by the GJP than by the MDP.
We present here a mathematical model that integrates various bystander responses, such as MDP and GJP signaling, DSB induction, cell-cycle arrest, and cell death. Because it simulates spatial and temporal conditions of irradiation and cellular characteristics, our model will be a powerful tool to predict dynamical radiobiological responses of a cellular population in which irradiated and non-irradiated cells co-exist.
It has long been desired to understand low-dose effects of ionizing radiation on living systems. From a practical viewpoint, the risks of low-dose radiation have been discussed for clinical diagnosis and therapy, as well as for radiation control around a nuclear disaster site. Even though only a limited number of cells in the population are exposed to doses of radiation below 1 Gy, it has been widely recognized that radiobiological effects, such as chromosome aberrations  and micronucleation or apoptosis , can be induced not only in the irradiated cells but also in non-irradiated cells. The transfer of radiation effects from irradiated to non-irradiated cells is termed the “bystander effect”. Bystander effects are activated by two intercellular signaling pathways: the medium-mediated pathway  and the gap junctional pathway . To understand the mechanisms of various bystander responses, the detailed dynamics of these two pathways is needed.
One powerful approach to understand the mechanism of bystander signal transfer is modeling the processes mathematically. In the last decade, several studies have reported models of the bystander effect. A Bystander and Direct (BaD) model, based on surviving fractions, was proposed by Brenner et al. . In the BaD model, cell death could be induced both by direct exposure to ionizing radiation and by bystander signals. Ebert et al. developed a surviving fraction model taking into account the density of bystander signals and the probabilities of interaction between a cell and the bystander signals . Several surviving fraction models or simulations have been reported that focus on the spatial and temporal kinetics of the bystander effect through the medium-mediated pathway [7, 8]. These models simulate the cellular population response as a surviving fraction and allow easy analysis of experimental data. However, it is very difficult to simulate the individual kinetics of cellular responses using these previous models. Khvostunov and Nikjoo developed a biophysical model that takes into account the diffusion of bystander signals and individual cell death (bystander diffusion modeling: BSDM) [9, 10]. In this model, the reactions of signals were calculated using the Monte Carlo method. Several similar attempts have also been made [11–13]. Richard et al. developed a simple deterministic model based on a cellular automaton [14, 15], which a useful model to describe the dynamics of a population of proliferating cells. Another model includes the repair of cell damage and cell death for both tumor cells and normal cells . Although these models could simulate individual cellular responses caused by signaling through the medium-mediated pathway, the gap junctional pathway has not been fully taken into account.
In some of these studies [5, 6, 9, 10, 12–15], biological end-point of bystander cells were simply described as binary states, namely, “life” or “death”. However, this is not enough to represent the complex responses of bystander cells, which could have various intermediate states between life and death, such as arrest or elongation of the cell cycle or the G0 state. It has also been widely recognized that cell-cycle arrest would be strongly dependent on the number of double-strand breaks (DSBs) in an irradiated cell, although there have been very few reports on bystander cell-cycle modifications .
We have previously developed a model framework for bystander effects on cell-cycle modification or cell death mediated by intercellular signaling through both the medium-mediated pathway and the gap junctional pathway . This model could calculate the time course of the individual cellular responses, however, cellular DNA damage was not fully taken into account, even cell-cycle arrest is thought to be one of the primary processes to secure sufficient time for repairing DNA damage, such as DNA DSBs.
In this study, we improved our model to simulate the induction of DSBs, cell-cycle modification, and cell death arising in non-irradiated cells. To model the inhomogeneous diffusions through the gap junctional pathway, we introduced “diffusion-direction constants” that define the direction of signal diffusion by referring to the states of adjacent cells. Furthermore, we have also introduced a biological “apparatus”, a virtual clock regulating cell cycles by cell-cycle checkpoints which function based on the number of DSBs produced. These are major advantages of our bystander model over previous ones. Using this model, we have analyzed in detail the dynamics of populations of cells, particularly cell-cycle modification or cell death. Various conditions, such as radiation dose or the relative contribution of each signaling pathway were carefully examined.
The simulation system is based on a two-dimensional cellular automaton, as shown in Fig. 1 a. The grids are categorized into three types: cell, medium, and wall. Each cell grid consists of a cell and culture medium. Each medium grid is filled entirely with culture medium. The wall grids represent a cell-culture container.
The simulation algorithm consists of four components: (1) irradiation of the cells, (2) generation of bystander signals in the irradiated cells and diffusion of these signals to surrounding cells, (3) cellular DNA damage induced by irradiation or bystander signals, and (4) cellular responses to the DNA damage. Cell-cycle modification and cell death were simulated as cellular responses. These four components relate to each other as shown in Fig. 1 b. The spatial and temporal evolution of the cellular population can be simulated by following the temporal progress of these components in individual cells.
Irradiation to cells
The biological effectiveness of ionizing irradiation generally depends on both absorbed dose and dose rate. The dose is usually described as the average dose calculated by dividing the total absorbed energy by the mass of the target cellular population. However, the radiation dose is really the result of the random traversal of radiation tracks and the absorbed dose in each cell is non-uniform, particularly following low-dose irradiation.
We therefore assumed that a single radiation track passing through a cell gives a fixed radiation dose, D 1track, which is set to correspond to an elementary dose. The number of radiation tracks in each cell is calculated based on random function, and the absorbed dose is calculated as the accumulation of the radiation tracks. The number of radiation tracks in grid position (i,j) at time t is represented by a random variable K i,j (t). Therefore, the absorbed dose, R i,j (t), in (i,j) at t is
The value of K i,j (t) is determined according to a Poisson distribution:
where K P i,j (K n) is the probability of K n radiation tracks arising in grid (i,j) over time interval Δ t and K a is the average number of radiation tracks passing through a grid in interval Δ t.
The values of D 1track and K a can be determined for various radiation types. For example, when cells are irradiated by 60Co γ-rays, D 1track=0.001 Gy. This value can be derived from the average dose given by secondary electrons generated by γ-irradiation at a very low dose rate . In the case of heavier ions, such as helium or carbon ions, D 1track can be higher.
Generation and diffusion of intercellular signaling
Intercellular signaling is described by the diffusion of transmitters into surrounding grids. We considered two kinds of radiation-induced signal: one that is transferred through the medium-mediated pathway (MDP) and another one through the gap junctional pathway (GJP). It has been mostly understood that some substances, such as NO radicals or cytokines, are major signals transferred through MDP [20, 21]. On the other hand, for the GJP, not all substances involved in the signaling have been identified. In general, cells communicate using messenger molecules such as c-AMP, calcium ions, inositol trisphosphate, ATP, or other nucleotides through gap-junctions [22–25] with a hexameric structure consisting of various connexin proteins . Among these chemicals, the calcium ion is known to play a major signaling role in various biological processes. Thus, we assume that the virtual signal through the GJP has the same characteristics as that of calcium ions in our model.
The two virtual signals were assumed to have the following properties: (1) there are no interactions between the signal molecules; (2) the signals are proportional to absorbed dose and degraded exponentially over time; and (3) the signals are diffusion-controlled. Previously reported bystander-signal molecules, such as NO radicals or cytokines, are not likely to specifically recognize and interact with each other like proteins, although they could be degraded non-specifically by radical reactions or the digestive action of cells. Therefore, we made the no interaction and exponential decay assumptions. Assumption (2) seemingly contradicts the experimental result reported by Hu et al. that there was no dose dependence of the number of cells possessing DSBs . In their study, bystander cells were counted as “DSB-positive cells” irrespective of the number of DSBs actually arising in the cell. Thus it is difficult to directly compare the number of bystander signals with the number of resulting DSBs. In order to avoid possibly underestimating the efficiency of the number of bystander signals or induced DSBs, assumption (2) is not unrealistic as a first order of approximation. Indeed, there has been further experimental evidence of the dose dependence of bystander cell death in the low-dose region . In higher-dose regions, however, there have been no reports of a dependence of produced signals on dose. To validate our simulation results, experimental investigations should be performed in future (see the ‘Discussion and conclusions’).
The diffusion equation is generally discretized in the following way (see Additional file 1):
where ϕ i,j (t) is the signal concentration in grid (i,j), Δ t is the time interval, d is the width of the grid, ϕ w is the diffusion coefficient, and (k,l) represents the positions of the eight grids surrounding grid (i,j). Based on the distance of (k,l) from (i,j), we categorized the surrounding grids into two parts, indexed by (k 1,l 1) and (k 2,l 2). The former contains grids located to the right, left, above and below the central grid, and the latter contains the skew grids. That is, (k 1,l 1) is (i+1,j), (i−1,j), (i,j+1), or (i,j−1), while (k 2,l 2) is (i+1,j+1), (i−1,j+1), (i+1,j+1), or (i−1,j−1).
Because the virtual signal diffuses into both cell grids and medium grids by the MDP, the signal is homogeneously transferred in all directions. On the other hand, the virtual signal through the GJP transferred to only nearby cell grids, resulting in an inhomogeneous diffusion of the signal (shown in Fig. 2). To control the direction of signal diffusion, we defined diffusion-direction constants, M w k,l′ and G w k,l′:
where M w and G w are diffusion constants. Here, we note that the cells are in a three dimensional condition of cultured dish. The volume of medium is much larger than the total volume of those of cells attached to the bottom of the dish, so the diffusion constant of the MDP in a cell grid was set to the same value as that for a medium grid. The diffusion-direction constants show the direction of intercellular signaling (red and blue arrows in Fig. 2). When the grid (k,l) is a cell grid, M w k,l′ and G w k,l′ are set as diffusion coefficients leading signals in (i,j) to grid (k,l). When the grid (k,l) is a medium grid, G w k,l′ is set to “0”, indicating that the signal through the GJP does not diffuse into medium grid (k,l). In contrast, the signal through the MDP diffuses into grid (k,l). There is, of course, also no diffusion into wall grids.
The quantities of the virtual signals through the MDP and the GJP are described by M i,j (t) and G i,j (t), respectively (red and blue circles in Fig. 2). Based on Eqs. (3–5), M i,j (t) and G i,j (t) in (i,j) at t are given by
where M α and G α are signal-production constants, and M α R i,j (t) and G α R i,j (t) are quantities of the produced virtual signals through the MDP and the GJP, respectively. R i,j (t) is the absorbed dose calculated based on the number of radiation tracks (Eq. (1)). The radiation track produces constant quantities of virtual signals represented by a “unit”. The parameters M β and G β are decay constants, and M β M i,j (t) and G β G i,j (t) are reductions in quantities of virtual signals, for the MDP and the GJP, respectively.
Production of DSBs
DNA DSBs are one type of severe damage induced by irradiation and causes cell-cycle arrest or cell death. Bystander signals are also known to induce DSBs. Several studies have shown that the number of DSBs is an important factor in cell-cycle modification. For example, when the number of DSBs is smaller than a certain threshold, cells are released from G2 arrest . To determine the cellular response to irradiation and also bystander signals, we estimated the number of DSBs. We made two assumptions to model DSB production: (1) the radiation and the bystander signals induce DSBs independently; (2) the number of produced DSBs is proportional to the quantities of bystander signals. It is known that the number of DSBs is proportional to the absorbed radiation dose . Because the relationship between the quantity of signal and DSB production has not been clarified, we assumed a linear relationship. Based on assumption (1), the number of DSBs, Z i,j (t), is represented by
where the number of DSBs induced by radiation, virtual signals through the MDP and the GJP are R Z i,j (t), M Z i,j (t) and G Z i,j (t), respectively. The number of DSBs endogenously induced by background factors such as reactive oxygen species is represented as B Z i,j (t)  and the number of repaired DSBs is represented as r Z i,j (t).
Since low-absorbed-dose irradiation induces a small number of DSBs , the distributions of R Z i,j (t), M Z i,j (t), G Z i,j (t), and B Z i,j (t) are based on Poisson distributions. The distribution of R Z i,j (t) is
where ZR P i,j (ZR n) is the probability of ZR n DSBs induced by radiation arising in a cell over an interval Δ t. ZR a i,j (t) is the average of R Z i,j (t), and is calculated as ZR λ i,j R i,j (t) because of the proportional relationship between the number of DSBs and the absorbed radiation dose . Here, ZR λ i,j is the induction coefficient for DSBs induced by irradiation.
Similarly the distributions of M Z i,j (t) and G Z i,j (t) are given by
where ZM P i,j (ZM n) is the probability of ZM n DSBs induced by the MDP arising in a cell over interval Δ t. ZG P i,j (ZG n) is the probability of ZG n DSBs induced by the GJP arising in a cell over Δ t. ZM a i,j (t) and ZG a i,j (t) are the averages of M Z i,j (t) and G Z i,j (t), and are calculated as ZM λ i,j M i,j (t)Δ t and ZG λ i,j G i,j (t)Δ t, respectively, based on assumption (2). Here, ZM λ i,j and ZG λ i,j are induction coefficients for DSBs induced by virtual signals through the MDP and the GJP, respectively.
The distribution of B Z i,j (t) is given by
where ZB P i,j (ZB n) is the probability of ZB n DSBs induced by background factors arising in a cell over Δ t. ZB a i,j is the average of B Z i,j (t) and ZB λ i,j is the corresponding induction coefficient.
The number of repaired DSBs, r Z i,j (t), is calculated using a probability of DSB repair. Most of the produced DSBs are known to be repaired in living cells through several independent pathways. Radiation-induced DSBs are also thought to have various complexities  and to be repaired with different kinetics in each pathway. Although some repair kinetics models have been reported in previous studies [34, 35], there is little experimental evidence that reaction rate constants are completely discriminated for each pathway. The relationship between DSB complexity and repair pathway has also not been elucidated. Hence, in the current stage of our model, we assumed for simplicity that all DSBs are repaired with the same probability, Zr λ i,j .
The calculation algorithm for r Z i,j (t) is shown in Fig. 3. First, r Z i,j (t) is set to 0. The index z in the algorithm (Fig. 3) counts the number of DSBs, and is initially set to 0. When Z i,j (t−Δ t) is larger than z, a uniform random number, r P, in [0,1], is generated. When r P> Zr λ i,j , the DSB is not repaired. When r P is smaller than Zr λ i,j , the DSB is repaired and r Z i,j (t) is increased by one. After the comparison between r P and Zr λ i,j , z is increased by one. The generation of r P and the comparison are repeated until z reaches Z i,j (t−Δ t).
The parameters ZR λ i,j , ZM λ i,j , ZG λ i,j , ZB λ i,j , and Zr λ i,j are initially set to different values for individual grids. To reflect the characteristics of individual cells, we assume that the parameters are taken from the positive part of a normal distribution.
Cell-cycle arrest is known to occur at certain checkpoints when DNA is damaged, and modification of the cell cycle is an important index to measure when monitoring radiation-induced responses. However, radiation-induced cellular responses have been estimated mainly based on cell death so far. In our model, we consider both cell cycle progression and cell death after irradiation.
The phase of the cell cycle or cell death for the cell grid (i,j) at time t is represented by S i,j (t). The progress or arrest of the cell cycle is determined by a virtual clock, as shown in Fig. 4. The virtual clock is divided into five phases, G1, S, G2, M1, and M2. M1, and M2 are a substructure of the M phase, which is divided at the M1/M2 checkpoint. The checkpoints (G1/S, S/G2, G2/M1, and M1/M2) are the time points defined by the clock. The time of cell division is represented by the M2/G1 time point. S i,j (t) is determined by the position of the virtual clock hand. When the clock hand stops at each checkpoint, S i,j (t) indicates cell-cycle arrest.
Cell death is generally divided into reproductive death  and interphase death . Reproductive death is the loss of the proliferative ability of the cell, and cells keep their cellular activity even after stopping cell division. Interphase death shows no proliferation, and the cells are disrupted. We modeled both types of cell death, taking into account that the reproductively dead cells still transfer signals through the GJP.
Cellular states are represented by four states, the proliferating (PR), pre-reproductive death (p-RD), reproductive death (RD), and pre-interphase death (p-ID) states, as shown in Fig. 5. Each state has a virtual clock. We used S i,j (t) to represent not only the phase of the cell cycle but also these four states. In the PR state, the cells have infinite proliferative capacity. In the p-RD state, the number of cell divisions is randomly assigned from one to three. When the cell division stops in the p-RD state, the cell transits to the RD state, in which the cell remains but the virtual clock stops permanently. In the p-ID state, the number of cell divisions is less than one and so the cell disappears. The transition pathways from the PR state to the p-RD or p-ID states are assumed to be irreversible.
It has been reported that the G2/M checkpoint has a threshold in the number of DSBs to determine the cell-cycle progression or arrest . In our model, we set threshold values using the number of DSBs, Z i,j (t). The threshold for transition between the arrest and the progression of a cell cycle was represented by A H i,j . The thresholds for transition from the PR to the p-RD state was represented by p−RD H i,j , and that for transition from the PR to the p-ID state by p−ID H i,j . The state of the cell, S i,j (t), is in progression or is arrested when Z i,j (t) lies below or exceeds the threshold values for each phase of the cell cycle. Figure 6 schematically shows the threshold values determining the state transitions. When Z i,j (t) is larger than the threshold, p−RD H i,j , for transition from the PR to the p-RD states, S i,j (t) takes the “p-RD:G1”. Similarly when Z i,j (t) is larger than the threshold, A H i,j , for transition to G1 arrest, S i,j (t) takes the “p-RD:G1 arrest” and the clock hand stops. When Z i,j (t) is smaller than A H i,j , Z i,j (t) is in the “p-RD:S” state and the clock hand is progressing.
When S i,j (t) is the PR or the p-RD state and reaches the M2/G1 checkpoint, the cell is divided into two daughter cells. One of the adjacent medium grids is replaced by one of the daughter cells. If there is no medium grid in the adjacent area, S i,j (t) takes the G0 phase.
The parameter values for state transitions A H i,j , p−RD H i,j , and p−ID H i,j are set differently for each individual grid. To reflect the characteristics of individual cells, we assumed that the parameters are taken from the positive part of a normal distribution. All the variable numbers and parameters used in our model are shown in Tables 1 and 2.
To simulate individual cellular responses by bystander signals, our model needs to set the values of the 19 parameters shown in Table 2. Some of these values can be found from previous studies. We estimated the others using experimental data as follows.
As mentioned above, some of parameters have a distribution which was generated by the algorithm in Fig. 7. First, a parameter value, V i,j , is generated using a function for producing a normal random number, RAND f(AV R, SD R), where AV R and SD R are the average and standard deviation of the normal random number. Next, V i,j is compared with the minimum value of the distribution, MIN V. Until the condition V i,j > MIN V is satisfied, V i,j is regenerated by RAND f(AV R, SD R). The parameter values of all grids were determined using this method. To estimate AV R and SD R for each parameter, we fitted the calculated data to selected experimental data. For the calculation of cell-cycle modification by bystander effects, there is no parameter set experimentally obtained from one cell line. Thus, we assembled parameter values from various cell lines used in previous studies.
Change in the number of DSBs with time
To estimate AV R and SD R for the probability of DNA repair, Zr λ i,j , data calculated by our model were fitted to a time course of the number of DSBs determined by the number of γ-H2AX foci per cell after X-ray irradiation . The experimental data showed that the number of DSBs after irradiation (1.5 Gy) decreased with time.
We calculated the number of DSBs using Eq. (8). So many DSBs were produced by irradiation in the experiment that the number of DSBs caused by background factors or bystander signals could be neglected. Thus, the number of DSBs induced by virtual signals through the MDP, M Z i,j (t), and through the GJP, G Z i,j (t), were both set to 0. The experimental data used for data fitting (shown in Fig. 8 a) indicated that the number of DSBs just after (5 min) irradiation (1.5 Gy) was 30–35 . We therefore set R Z i,j (t=0)=35. Because the number of DSBs endogenously induced by background factors, B Z i,j (t), was much smaller than R Z i,j (0), B Z i,j (t) was set to 0. To calculate the average number of DSBs, we used 10,000 grids. We performed a preliminary search and found the ranges of parameter values that could reproduce the experimental data. The ranges for AV R and SD R were set to be from 1.0×10−8 to 1.0×10−6 and from 0 to 1.0×10−4, respectively. We determined the values of AV R and SD R of Zr λ i,j in the G1 phase by data fitting to be 9.33×10−7 s −1 and 2.0×10−5 s −1, respectively, and the results are shown in Fig. 8 a.
Number of cells that have DSBs following intercellular signaling
To estimate AV R and SD R for the induction coefficients for DSBs induced by virtual signals and background factors, ZM λ i,j , ZG λ i,j , and ZB λ i,j , results calculated using our model were fitted to the number of bystander cells with DSBs after α-irradiation . Since the irradiated dishes in the experiment were filled with cells in the G0 phase, we used a space composed of 100 ×100 grids in the G0 phase. The grid width, d, was set to 10 μm.
First, we estimated AV R for ZB λ i,j using control data in the experiment . The SD R of ZB λ i,j was set to 0 since the standard deviation is much smaller than the average. Because the control data have been measured under sham-irradiation conditions, the absorbed radiation dose and the quantities of virtual signals through the MDP and the GJP, R i,j (t), M i,j (t), and G i,j (t), were set to 0. We then used Eqs. (8) and (15) to calculate the number of DSBs, Z i,j (t). The probability of DNA repair, Zr λ i,j , was set to the value for the G0 phase estimated by the method described in the previous section. In the previous report, cells with one or more γ-H2AX foci were scored as cells with DSBs . Based on this experimental method, we calculated Z i,j (t) in 10,000 grids over 30 min, then the number of cells that had Z i,j (t)>0 were counted as cells with DSBs. The range of searching for the value of AV R for ZB λ i,j was set to be from 5.0×10−6 s −1 to 1.5×10−5 s −1. We determined values of AV R for ZB λ i,j in the G0 phase by data fitting, and the results are shown in Fig. 9. The value of AV R was estimated to be 1.4×10−5 s −1. The value of AV R for ZB λ i,j in the G1 phase was set to be the same as in the G0 phase. The values of AV R for ZB λ i,j in the S, G2, M1, and M2 phases were set to 2.8×10−5 s −1 because the amount of DNA in these phases is double that in the G1 phase.
For ZM λ i,j , we used the data  of cells treated with lindane which was used as an inhibitor of gap-junctional intercellular signaling, and the virtual signal through the GJP, G i,j (t), as well as the absorbed dose, R i,j (t), were set to 0. Then we used Eqs. (6), (8), (11), (12) and (15) to calculate the number of DSBs, Z i,j (t). We assumed that the virtual signal through the MDP had the same characteristics as those for cytokines, such as interleukin 8. The signal-production constant of the virtual signal through the MDP, M α, was set to 1 unit/Gy where “unit” denotes an arbitrary quantity of the virtual signal. Because the bystander signals through the MDP were still active at 60 hours after irradiation , the decay constant of the virtual signal through MDP, M β, was set to 4.6×10−6 s −1. The diffusion constant of the virtual signal through the MDP was set to 1.0×10−10 m 2/s . The initial value of virtual signals through the MDP, M i,j (t), in all cell grids was set to 2.5×10−3 units based on the radiation dose and the number of irradiated cells in the experiment. The values of ZB λ i,j and Zr λ i,j were set to the estimates shown in Table 2. The ranges of searching for the values of AV R and SD R were set to be from 2.0×10−3 to 4.0×10−3 DSBs/unit/s and from 0 to 1.0×10−4 DSBs/unit/s, respectively. We determined the values of AV R and SD R by data fitting, and the results are shown in Fig. 9. The values of AV R and SD R were estimated to be 6.0×10−3 DSBs/unit/s and 1.0×10−2 DSBs/unit/s, respectively.
For ZG λ i,j , we used data obtained at 30 min after 1 cGy-irradiation  and used Eqs. (6–8) and (11–15) to calculate the number of DSBs, Z i,j (t). As mentioned above, the virtual signal through the GJP was assumed to have the characteristics of calcium ions. The signal-production constant of the virtual signal through the GJP, G α, was set to 1 unit/Gy. The decay constant, G β, was set based on the efflux of calcium ions from the cell. The velocity of the efflux of calcium ions could be divided into a fast phase and a slow phase . However, the life time of the fast phase signal was too short to affect cellular responses ranging from minutes to hours, so we ruled out this phase. The value of G β was set to 1.18×10−3 s−1 based on the rate constant of the slow phase. The diffusion constant of the virtual signal through the GJP was set to 5.0×10−11 m 2/s, consistent with reported values . The values of ZM λ i,j , ZB λ i,j , and Zr λ i,j were set to the estimates shown in Table 2. The initial values of M i,j (t) and G i,j (t) were set to 2.5×10−3 units. The ranges for searching for the values of AV R and SD R were from 2.0×10−3 to 8.0×10−2 DSBs/unit/s and from 0 to 1.0×10−4 DSBs/unit/s, respectively. We determined the values of AV R and SD R by data fitting, and the results are shown in Fig. 9. The values AV R and the SD R were finally estimated to be 6.0×10−2 DSBs/unit/s and 1.0×10−2 DSBs/unit/s, respectively.
Cell survival fractions
The values of AV R and SD R for the induction coefficient for DSBs induced by radiation, ZR λ i,j , and the threshold of the DSB number for transition from the PR to the p-RD state, p−RD H i,j , were estimated based on cell survival fraction data. A cell grid in the PR state was defined to be equivalent to a surviving cell observed in experiments.
The survival of human cancer cells (HeLa cells) is the only data to include survival fractions for each cell-cycle phase . The parameter setting and the number of grids were determined in the same way as in the previous subsection. We used Eqs. (1–2) and (8–10) to calculate the number of DSBs, Z i,j (t). The expected value of the absorbed dose given by one radiation track, D 1track, was set to 1.0×10−3 Gy. The dose rate was set to 0.5 Gy/min. The irradiation time was set to 1, 2, 3, 4, 5, 6, 7, or 8 min according to the radiation doses in the experiment. The values of AV R for ZR λ i,j in the G0 and G1 phases were set to 40 DSBs/Gy because a previous study estimated that 1 Gy irradiation induced 40 DSBs . Since the amount of DNA in the S, G2, and M phases is double that in the G1 phase, the values of AV R in the S, G2, M1, and M2 phases were set to 80 DSBs/Gy. During the irradiation period, the number of induced DSBs was much larger than that of repaired DSBs. Thus, r Z i,j (t) was set to 0. We calculated the number of DSBs, Z i,j (t), in 10,000 grids. The number of surviving cells was defined as the number of grids that had Z i,j (t)< p−RD H i,j . In a similar way to the parameter estimation in previous sections, we determined the value of SD R for ZR λ i,j , and the values of AV R and SD R for p−RD H i,j by data fitting. The results are shown in Fig. 10. The determined values of AV R and SD R are shown in Table 2.
To analyze the dynamics and general properties of our model, we examined cell-cycle modification caused by intercellular signaling after exposure to various doses of radiation.
We used a simulation space composed of 103 ×303 grids, as shown in Fig. 11. The space was divided into three areas, A, B, and C, and cellular populations were introduced in each area. The cells in area A were cultured sparsely: cells occupied 15 % of a circle with a radius of 40 grids. In areas B and C, the cell populations were almost at 100 % confluence. Since the cells in the area A had enough space to progress their cell cycle for growth, the effect on cell-cycle modification is expected to be more easily observed in area A than in the other areas. In area B, only the central area of the population was irradiated. The non-irradiated cells received intercellular signals through both the MDP and the GJP, on the other hand the cells in areas A and C received signals through the MDP only. The contribution of the GJP to cell-cycle modification could be examined by comparing the populations in areas B and C.
Within the irradiation area, cells were irradiated uniformly (D 1track=1.0×10−3 Gy). The dose rate was set at 1 Gy/min and the tested doses were 0, 0.1, 0.2, 0.5, 1, 2, or 5 Gy. The case without exposure (0 Gy) was the control.
For the virtual signals through the MDP and the GJP, the parameters shown in Table 2 were used.
The number of DSBs that induced G2/M checkpoint release has been suggested to be approximately 20 . Thus, the average threshold, A H i,j , at G2/M1 checkpoints was set to 20 DSBs. The average thresholds were set to 5 and 40 DSBs for G1/S and S/G2, respectively (Shibata A., private communication, 2015). The M1/M2 checkpoint is the spindle checkpoint and does not directly relate to DSBs. Thus, the threshold value was set to 40, which was the largest value used in the present study. Standard deviations of the thresholds were set to 20 % of the averages.
For the averages and standard deviations of thresholds for transition from the PR to the p-RD state, p−RD H i,j , we used the values estimated by data fitting shown in Table 2.
The threshold of interphase death was determined to realize a high frequency of death in a higher dose region and low frequency in a lower dose region, as reported by Puck and Marcus . We assumed that the threshold, p−ID H i,j , was two times larger than the threshold for reproductive death based on the experimental frequencies of the two types of cell death.
An average cell cycle period was set to 24 h based on HeLa cell data . The periods of the G1, S, G2, M1, and M2 phases were set to 11, 8, 4, 0.5, and 0.5 h, respectively. The calculated cell culture time was set to 96 h to simulate cell fates even when the cycles were greatly prolonged.
We performed 10 trials using various initial conditions for the positions of cell grids and random numbers for the calculation of Poisson distributions.
The quantities of virtual signals in grids (i=1, 2,…,303,j=52) for 0.1 Gy irradiation are shown in Fig. 12, indicating the quantity of signal through the MDP, M i,j (t), (Fig. 12 a), and through the GJP, G i,j (t) (Fig. 12 b). At t=0, quantities of virtual signals at the irradiation area were around 1 unit. After 5 min, M i,52(t) for the population in area B was smaller than G i,52(t), reflecting the size of the diffusion area. After 10–60 min, the virtual signals through the MDP were transmitted to distant grids and finally diffused into the whole area. On the other hand those through the GJP were transmitted within area B, and finally almost disappeared. These results indicate that non-irradiated cells could be exposed to an MDP signal over a longer period than a GJP signal.
Dependence of cell-cycle modification on radiation dose
We examined the dynamics of the population in all the areas, as shown in Figs. 13 and 14. The irradiated cells were excluded from the cell counting in area B, and the cells were subtracted from the number of cells at G0 in area C to compare these two areas. Time courses for 0.1, 0.2, and 0.5 Gy irradiation were similar to that of the control. However, those for higher doses (1, 2, and 5 Gy) significantly decreased with increasing dose. The numbers of cells in each cell-cycle arrest mode after irradiation is shown in Figs. 15 and 16. In the low-dose region (0.1–0.2 Gy) there were no significant effects for all cases. However, above 0.5 Gy, the number of arrested cells increased with dose in all cases except for G1 arrest. Area A showed significant cell-cycle modification by intercellular signaling for all cases except for G0.
Area A was analyzed further to investigate dose dependence. The number cells decreased with time for the G1, S, and G2 phases, but increased for the G0 phase. For 1 Gy irradiation, the increases of cells in all phases were delayed relative to the control. For 5 Gy irradiation, cells in all phases except G0 significantly decreased with time, and the number of G0 cells remained extremely low throughout the simulation period. The number of G1 arrested cells after 0.1–2 Gy irradiation slightly increased, compared with that for the control. Due to DSB repair, the cells arrested in the G1 phase by exposure to 0.1–2 Gy irradiation were released from the G1/S checkpoint and entered back into the cell cycle. Thus the number of G1 arrested cells decreased with time. For exposure to 5 Gy irradiation, the number of arrested cells increased greatly with different kinetics for each arrest mode: the rate of increase strongly correlated with the threshold values. Because almost all cells in the case of 5 Gy irradiation showed cell-cycle arrest and the number of G1 cells decreased, the number of G1 arrested cells was smaller than that in the case of 2 Gy irradiation. The time courses of population growth for each radiation dose are shown in Fig. 17. In area A, the population growth for 0.1–0.5 Gy exposure was similar to that of the control. For the higher dose range of 1–5 Gy, the rate of increase was suppressed because almost all cells were arrested by the intercellular signaling. From Fig. 17 b, it can be seen that the number of cells decreased with dose but the decrease was not proportional to dose.
Contribution of each intercellular signaling pathway
We now focus on areas B and C after irradiation (Figs. 13, 14, 15, 16 and 17). We scored cells with the exception of irradiated cells. In all cell-cycle phases, the numbers of cells in both areas decreased with time when exposed to 5 Gy (Figs. 13 and 14). In particular the number of cells strikingly decreased in the area B compared to area C due to the quantities of virtual signals in the two pathways. The cells in area B received virtual signals through both the MDP and the GJP, while the cells in area C received signals only through the MDP.
From Figs. 15 and 16, it can be seen that the total numbers of arrested cells for the control were negligible. For 0.5–5 Gy exposure, the number of G1 arrested cells in area B increased with time, and was similar to that in area C except for 5 Gy. The numbers of S, G2, and M arrested cells in area B for 1–5 Gy doses was smaller than those in area C: arrested cells were more efficiently accumulated in the population of area B larger than in the population of area C. From Fig. 17, it can be seen that the time course and the dose response of populations in area B were similar to those of populations in area C.
Discussion and conclusions
To investigate the dynamics of a cell population, the numbers of cells in the PR state, the pre-reproductive death (p-RD) state, the reproductive death (RD) state, and the pre-interphase death (p-ID) state were analyzed as shown in Fig. 18. The number of p-RD state cells reached a maximum 60–70 h after irradiation because the cells transited to the RD state after 24–72 h (1 to 3 cell divisions). The number of p-ID state cells reached a maximum 20–30 h after irradiation because the cells transited to interphase death within 1 cell division. The number of PR state cells decreased with dose because of cell-cycle arrest or cell death. During the 0–50 h period, the number of RD state cells was much smaller than the numbers of G1, S, and G2 arrested cells (Figs. 13 and 14). Thus, the cell-cycle arrest induced by bystander signals reduced the number of PR state cells soon after irradiation.
The number of arrested cells did not increase linearly with dose. This is due to the DSB thresholds determining cell-cycle arrest. Since we set the threshold at the G1/S checkpoint to be smaller than that at other checkpoints, G1 arrested cells accumulated in the population more than other arrested cells. Although the threshold at the G2/M checkpoint was set to half of that at the S/G2 checkpoint, the number of G2 arrested cells was more than twice the number of S arrested cells. Thus, the number of arrested cells depends sensitively on these thresholds.
To validate the simulation results, the dynamics of the model should be further examined by comparison with applicable experimental data. Although there have been few reports on bystander cell-cycle effects so far because of technical difficulties in tracking the cell cycle of individual cells, recent live cell imaging studies using fluorescent ubiquitination-based cell-cycle indicators (FUCCI) have enabled visualization of the cell cycle of individual cells . Such an approach is expected soon to provide applicable data that can be used to refine and improve our model.
Using FUCCI cells, we recently determined cell-cycle periods after exposure of specific cells in a population to X-rays  or an X-ray microbeam . To validate the present simulation results, particularly Figs. 13, 14, 15 and 16, we have also been examining whether non-irradiated cells undergo cell cycles in a colony in which some cells were exposed to an X-ray microbeam (to be submitted). Thus, the present model could be used to predict dynamic changes from cycle-progressing cells to arrested cells, or the growth of cellular populations before performing irradiation experiments.
The GJP contributes more to bystander effects than the MDP, as shown in Figs. 13, 14, 15 and 16. For 5 Gy irradiation, the average number of DSBs in area B reached a maximum at 1 h after irradiation, and the maximum was 47 DSBs/cell in the G1 phase. In area C, the average value reached a maximum at 10 h, and the maximum was 17 DSBs/cell in the G1 phase. These results clearly show that the GJP contributed significantly to the induction of DSBs. Furthermore, the arrested cells were more efficiently accumulated in the population of area B than in the population of area C, indicating that signaling through the GJP plays a significant role in bystander cell-cycle effects.
Experimentally suppressing all MDP or GJP signals would highlight the particular contribution of a specific pathway; however, it would be technically difficult to completely quench the signals using high concentrations of scavenging chemicals, which might cause undesirable side-effects in the cells. The present model is therefore a powerful tool for predicting bystander effects from the actions of intercellular signaling pathways.
Gap junctional pathway (intercellular signaling pathway)
Medium-mediated pathway (intercellular signaling pathway)
Nagasawa H, Little JB. Induction of sister chromatid exchanges by extremely low doses of alpha particles. Cancer Res. 1992; 52(22):6394–6.
Prise KM, Belyakov OV, Folkard M, Michael BD. Studies of bystander effects in human fibroblasts using a charged particle microbeam. Int J Radiat Biol. 1998; 74(6):793–8.
Mothersill C, Seymour C. Medium from irradiated human epithelial cells but not human fibroblasts reduces the clonogenic survival of unirradiated cells. Int J Radiat Biol. 1997; 71(4):421–7.
Azzam EI, de Toledo SM, Little JB. Direct evidence for the participation of gap-junction mediated intercellular communication in the transmission of damage signals from alpha-particle irradiated to nonirradiated cells. Proc Natl Acad Sci U S A. 2001; 98(2):473–8.
Brenner DJ, Little JB, Sachs RK. The bystander effect in radiation oncogenesis: II. A quantitative model. Radiat Res. 2001; 155(3):402–8.
Ebert MA, Suchowerska N, Jackson MA, McKenzie DR. A mathematical framework for separating the direct and bystander components of cellular radiation response. Acta Oncol. 2010; 49(8):1334–43.
Shuryak I, Sachs RK, Brenner DJ. Biophysical models of radiation bystander effects: 1. Spatial effects in three dimensional tissues. Radiat Res. 2007; 168(6):741–9.
McMahon SJ, Butterworth KT, Trainor C, McGarry CK, O’Sullivan JM, Schettino G, et al. A kinetic-based model of radiation induced intercellular signalling. PLOS One. 2013; 8(1):e54526.
Khvostunov IK, Nikjoo H. Computer modelling of radiation-induced bystander effect. J Radiol Prot. 2002; 22(3A):A33–7.
Nikjoo H, Khvostunov IK. Biophysical model of the radiation-induced bystander effect. Int J Radiat Biol. 2003; 79(1):43–52.
Ballarini F, Alloni D, Facoetti A, Mairani A, Nano R, Ottolenghi A. Modelling radiation-induced bystander effect and cellular communication. Radiat Prot Dosimetry. 2006; 122(1–4):244–51.
Xia J, Liu L, Xue J, Wang Y, Wu L. Modeling of radiation-induced bystander effect using Monte Carlo methods. Nucl Instrum Methods Phys Res B. 2009; 267(6):1015–8.
Sasaki K, Wakui K, Tsutsumi K, Itoh A, Date H. A simulation study of the radiation-induced bystander effect Modeling with stochastically defined signal reemission. Comp Math Methods Med. 2012; 2012(55):1–5. doi:10.1155/2012/389095.
Richard M, Webb RP, Kirkby KJ, Kirkby NF. computer model of the bystander effect: Effects of individual behaviours on the population response. Appl Radiat Isot. 2009; 67(3):440–2.
Richard M, Kirkby KJ, Webb RP, Kirkby NF. Cellular automaton model of cell response to targeted radiation. Appl Radiat Isot. 2009; 67(3):443–6.
Powathil GG, Munro AJ, Chaplain MAJ, Swat M. Bystander effects and their implications for clinical radiation therapy Insights from multiscale in silico experiments. Quant Methods. 2014. arXiv:1407.0867.
Kaminaga K, Noguchi M, Narita A, Sakamoto Y, Kanari Y, Yokoya A. Visualization of cell cycle modification by X-irradiation of single HeLa cells using fluorescent ubiquitination-based cell cycle indicators. Radiat Prot Dosim. 2015; 166(1–4):91–4. doi:10.1093/rpd/ncv168.
Hattori Y, SUzuki M, Funayama T, Kobayashi Y, Yokoya A, Watanabe R. A mathematical modelof radiation-induced responses in acellular population including cell-to-cell communications. Radiat Prot Dosim. 2015; 166(1–4):142–7. doi:10.1093/rpd/ncv149.
Booz J, Feinendegen LE. A microdosimetric understanding of lowdose radiation effects. Int J Radiat Biol Relat Stud Phys Chem Med. 1988; 53(1):13–21.
Hei TK, Zhou H, Chai Y, Ponnaiya B, Ivanov VN. Radiation induced non-targeted response: mechanism and potential clinical implications. Curr Mol Pharmacol. 2011; 4(2):96–105.
Facoetti A, Ballarini F, Cherubini R, Gerardi S, Nano R, Ottolenghi A, et al. Gamma ray-induced bystander effect in tumour glioblastoma cells: a specific study on cell survival, cytokine release and cytokine receptors. Radiat Prot Dosimetry. 2006; 122(1–4):271–4.
Lawrence TS, Beers WH, Gilula NB. Transmission of hormonal stimulation by cell-to-cell communication. Nature. 1978; 272(5653):501–6.
Loewenstein WR. Junctional intercellular communication: the cell-to-cell membrane channel. Physiol Rev. 1981; 61(4):829–913.
Sáez JC1, Connor JA, Spray DC, Bennett MV. Hepatocyte gap junctions are permeable to the second messenger, inositol 1, 4, 5-trisphosphate, and to calcium ions. Proc Natl Acad Sci U S A. 1989; 86(8):2708–12.
Pearson RA, Dale N, Llaudet E, Mobbs P. ATP released via gap junction hemichannels from the pigment epithelium regulates neural retinal progenitor proliferation. Neuron. 2005; 46(5):731–44.
Kumar NM, Gilula NB. Molecular biology and genetics of gap junction channels. Semin Cell Biol. 1992; 3(1):3–16.
Hu B, Wu L, Han W, Zhang L, Chen S, Xu A, et al. The time and spatial effects of bystander response in mammalian cells induced by low dose radiation. Carcinogenesis. 2006; 27(2):245–51.
Schettino G, Folkard M, Prise KM, Vojnovic B, Held KD, Michael BD. Low-dose studies of bystander cell killing with targeted soft X rays. Radiat Res. 2003; 160(5):505–11.
Deckbar D, Birraux J, Krempler A, Tchouandong L, Beucher A, Walker S, et al. Chromosome breakage after G2 checkpoint release. J Cell Biol. 2007; 76(6):749–755.
Ward JF. Some biochemical consequences of the spatial distribution of ionizing radiation-produced free radicals. Radiat Res. 1981; 86(2):185–95.
Swenberg JA, Lu K, Moeller BC, Gao L, Upton PB, Nakamura J, et al. Endogenous versus exogenous DNA adducts: their role in carcinogenesis, epidemiology, and risk assessment. Toxicol Sci. 2011; 120(Suppl 1):S130–45.
Ojima M, Ban N, Kai M. DNA double-strand breaks induced by very low X-ray doses are largely due to bystander effects. Radiat Res. 2008; 170(3):365–71.
Nikjoo H, O’Neill P, Terrissol M, Goodhead DT. Computational approach for determining the spectrum of DNA damage induced by ionizing radiation. Radiat Res. 2001; 156:577–83.
Taleei R, Nikjoo H. Biochemical DSB-repair model for mammalian cells in G1 and early S phases of the cell cycle Repair of the double-strand breaks induced by low energy electrons: a modelling approach. Mutat Res. 2013; 756(1–2):206–12.
Friedland W, Kundrat P, Jacob P. Stochastic modelling of DSB repair after photon and ion irradiation. Int J Radiat Biol. 2012; 88(1–2):129–36.
Lea DE. Actions of Radiations on Living Cells, 2nd ed. London: Cambridge University Press; 1956.
Kelly LS. Radiosensitivity of biochemical processes. Brookhaven Symp Biol. 1961; 14:32–52.
Jacobson K, Wojcieszyn J. The translational mobility of substances within the cytoplasmic matrix. Proc Natl Acad Sci USA. 1984; 81(21):6747–51.
Borle AB. Kinetic analyses of calcium movements in HeLa cell cultures II. Calcium efflux. J Gen Physiol. 1969; 53(1):57–69.
Bathany C, Beahm D, Felske JD, Sachs F, Hua SZ. A high throughput assay of diffusion through Cx43 gap junction channels with a microfluidic chip. Anal Chem. 2011; 83(3):933–939.
Terashima T, Tolmach LJ. Variations in several responses of HeLa cells to X-irradiation during the division cycle. Biophys J. 1963; 3(1):11–33.
Goodhead DT. Initial events in the cellular effects of ionizing radiations: clustered damage in DNA. Int J Radiat Biol. 1994; 65(1):7–17.
Puck TT, Marcus PI. Action of x-rays on mamMalian cells. J Exp Med. 1956; 103(5):653–66.
Hall EJ. Radiobiology for the Radiologist. LIPPINCOTT WILLIAMS & WILKINS, a WOLTERS KLUWER business, Two Commerce Square, 2001 Market Street, Philadelphia, PA 19103 USA; 2012.
Sakaue-Sawano A, Kurokawa H, Morimura T, Hanyu A, Hama H, Osawa H, et al. Visualizing spatiotemporal dynamics of multicellular cell-cycle progression. Cell. 2008; 132(3):487–98.
Narita A, Kaminaga K, Yokoya A, Noguchi M, Kobayashi K, Usami N, et al. Real-time observation of irradiated HeLa-cell modified by fluorescent ubiquitination-based cell-cycle indicator using synchrotron X-ray microbeam. Radiat Prot Dosim. 2015; 166(1–4):192–196. doi:10.1093/rpd/ncv156.
We wish to acknowledge valuable discussions with Dr. Keiji Suzuki (Nagasaki University). The authors are also grateful to Dr. Atsushi Shibata (Gunma University) for comments on this paper. This work was supported by JSPS KAKENHI Grant Number 15K16088.
The authors declare that they have no competing interests.
YH conceived and designed the research, analyzed the data, developed the mathematical model, and wrote the manuscript. AY designed research and wrote the manuscript. RW designed research and wrote the manuscript. All authors read and approved the final manuscript.
This file contains the method for discretization of the diffusion equation and Figure. Modeling of the intercellular signaling was based on the discretized diffusion equation. (PDF 65.9 kb)
About this article
Cite this article
Hattori, Y., Yokoya, A. & Watanabe, R. Cellular automaton-based model for radiation-induced bystander effects. BMC Syst Biol 9, 90 (2015). https://0-doi-org.brum.beds.ac.uk/10.1186/s12918-015-0235-2
- Cell cycle
- DNA double strand break
- Intercellular signaling
- Radiation-induced bystander effect