- Email: [email protected]
Cellular automaton model of cell response to targeted radiation
ARTICLE IN PRESS Applied Radiation and Isotopes 67 (2009) 443–446
Contents lists available at ScienceDirect
Applied Radiation and Isotopes journal homepage: www.elsevier.com/locate/apradiso
Cellular automaton model of cell response to targeted radiation M. Richard a,b, K.J. Kirkby a, R.P. Webb a, N.F. Kirkby b, a b
Surrey Ion Beam Centre, Advanced Technology Institute, University of Surrey, Guildford, GU2 7XH, UK Centre for Fluids and Systems, School of Engineering, University of Surrey, Guildford, GU2 7XH, UK
abstract It has been shown that the response of cells to low doses of radiation is not linear and cannot be accurately extrapolated from the high dose response. To investigate possible mechanisms involved in the behaviour of cells under very low doses of radiation, a cellular automaton (CA) model was created. The diffusion and consumption of glucose in the culture dish were computed in parallel to the growth of cells. A new model for calculating survival probability was introduced; the communication between targeted and non-targeted cells was also included. Early results on the response of non-confluent cells to targeted irradiation showed the capability of the model to take account for the non-linear response in the low-dose domain. & 2008 Published by Elsevier Ltd.
1. Introduction
2. Model
The accurate measurement of survival of mammalian cells at low doses of radiation in the past 20 years has shown that it could not be extrapolated from the high dose response. Some cell lines show a hyper radio-sensitivity (HRS) to low doses (below 0.5 Gy), followed by an increased radio-resistance between 0.5 and 1 Gy (Joiner et al., 2001); this leads to survival fraction curves that cannot be described by the traditional linear-quadratic model. On the other hand, not only irradiated cells are affected by radiation, but also bystander non-irradiated cells in the neighbourhood; this has been termed as the bystander effect (BE) and was widely investigated using microbeams (Prise et al., 2005; Morgan, 2003). The mechanisms underlying the HRS and the BE are not fully understood yet. Different experimental procedures have been used and different conclusions, sometimes contradictory, have been drawn. Some models have been created to take into account either the HRS (Wouters and Skarsgard, 1994) or the BE (Little et al., 2005; Nikjoo and Khvostunov, 2003; Brenner et al., 2001). However, none of those models attempted to embed both phenomena. The work presented in this paper and in Richard et al. (submitted) includes the HRS and the BE, using the cellular automaton (CA) computing approach; the diffusion of glucose and the cell cycle are fully described. The response of V79 cells to targeted radiation is simulated.
A CA is a lattice of sites whose time evolution depends on its previous states and on interacting rules between a site and its neighbours. In this model, a first CA computed the diffusion of glucose and a second CA the cell population growth. A detailed description of the diffusion procedure and cell cycle model can be found in Richard et al. (submitted). 2.1. The diffusion of glucose The lattice represented the cellular environment (i.e. culture dish) in which the glucose was diffusing. The diffusion procedure, adapted from Bandini et al. (2001), consisted of two steps. Any site of the CA was divided into four (north, south, east and west) and glucose first diffused between two adjacent quarters of two neighbouring sites (see Fig. 1), following the equation xðt nþ1 Þ ¼ xðt n Þ ðxðt n Þ yðt n ÞÞ k;
D¼
Corresponding author. Tel.: +44 1483686577; fax: +44 1483686581.
E-mail address: [email protected] (N.F. Kirkby). 0969-8043/$ - see front matter & 2008 Published by Elsevier Ltd. doi:10.1016/j.apradiso.2008.06.044
kp0:5
(1)
In Eq. (1), x is the concentration of glucose in one area of a given site (e.g. south) and y the concentration of glucose in the adjacent area of the neighbouring site (e.g. north area of the south neighbouring site).The concentration of glucose is then homogenised on a single site of the lattice. If Dx is the space increment of the CA and Dt the time increment, and k the constant in Eq. (1), the diffusivity D of the glucose is given by
Dx2 k Dt
(2)
Growing cells consumed the diffusing glucose, so its concentration on the lattice decreased. The diffusivity of glucose was chosen to be 7.2 1010 m2/s on the surfaces of the dish free from cells
ARTICLE IN PRESS 444
M. Richard et al. / Applied Radiation and Isotopes 67 (2009) 443–446
(Beuling et al., 2000), and 3 1011 m2/s on the surfaces covered by cells (Casciari et al., 1988).
3
4 x(tn)=1
3
y(tn)=3
2
2
2
5 x(tn+1)=2 y(tn+1)=2
3 1
2.2. The cell cycle
4
The sites of the second CA represented biological cells and the CA evolution showed the growth of the population. The cell cycle was divided into mitosis (M phase) and interphase, i.e. first gap (G1 phase), synthesis (S phase) and second gap (G2 phase) (see Fig. 2). The lengths of M, S and G2 were fixed, while G1 had a probabilistic duration. A G1 cell was allowed to enter the S phase only if the amount of glucose it had absorbed since last division was greater than a given minimum (see Kirkby et al., 2002; Faraday et al., 2001). The consumption of glucose was computed by a Monod-like equation:
5 4
r¼
2 2
2
5 4
1
(3)
When the surrounding density of cells became too high or the surrounding concentration of glucose too low, the cells enter a quiescent state called G]. Values of the different parameters of the model are displayed in Table 1.
5
3
V m CðtÞ K m þ CðtÞ
Fig. 1. Diffusion process. Any site of the CA is divided into four regions, north, south, east and west adjacent to the corresponding neighbouring sites. The glucose first diffuses between the regions of a site and the adjacent regions of the neighbouring sites. The equation of diffusion is x(tn+1) ¼ x(tn)(x(tn))(y(tn)) k with k ¼ 0.5. The glucose concentration is averaged between the four regions of each site.
2.3. Time and space scales of the cellular automata The time scale of the glucose diffusion (order of second) was much smaller than the time scale of the evolution of the cell population (order of hour). As a consequence, it was not possible to choose the same time and space increment for both cellular automata. Therefore, the following relationship was decided:
DxDCA ¼ 5 DxBCA
In Eq. (4), DxDCA is the space increment of the diffusion CA (DCA) and DxBCA the space increment of the biological CA (BCA), which is equal to 12 mm. The time increment of the DCA was then deduced from Eq. (2) to be 0.625 s and the time increment of the BCA was taken to be 1 h. At any iteration of the program, the probabilities of updating either the BCA or the DCA were
G# if not enough nutrient or high cell density Radiation induced cell cycle arrest
M phase G1 phase
Radiation induced cell cycle arrest
G2 phase
pBCA Dt DCA ¼ pDCA Dt BCA
Check nutrient absorption before S Death if absorption
S phase
(4)
Fig. 2. The cell cycle. The cell cycle is divided into mitosis (M phase), during which cells divide, and interphase, during which cells prepare themselves for the next division. The interphase consists of gap 1 (G1 phase), synthesis of DNA (S phase) and gap 2 (G2 phase). The M phase is short compared to the interphase. Cells in G1 continuously check the amount of glucose absorbed since last division; they can enter S only if this amount is greater than an internal threshold. If cells fail to check this criterion after a maximum period of time (maximum length of G1), they die. Cells in G1 can also enter a non-cycling phase, G], if the cell density is too high or if the concentration of glucose is too low. Irradiated cells are arrested in G2, for DNA integrity checking.
(5)
In Eq. (5), pBCA and pDCA are the probabilities of updating the BCA and the DCA, respectively, and DtBCA and DtDCA are the time increments of the BCA and DCA, respectively. 2.4. Microbeam irradiation It was assumed that a cell was very sensitive to a low dose of radiation because the repair processes were not triggered; as the dose increased, those processes would be triggered and the cell would become less sensitive. The limit between the low- and high-dose domain, namely, the dose dc, was supposed to be
Table 1 Main parameters used in the two CA
Glioma V79
DxDCA (mm)
DxBCA (mm)
DtDCA (s)
DtBCA (h)
k
Vm (104 nmol/cell/h)
Km (104 nmol/site)
C0 (104 nmol/site)
60 60
12 12
0.625 0.625
1 1
0.5 0.5
4 40 (G1 phase) 4 (other phases)
75 75
400 400
DxDCA and DxBCA (resp. DtDCA and DtBCA) refer to the space (resp. time) increment of the diffusion and biological CA. k is the constant in Eq. (1) in the main text and C0 is the initial concentration (in mol per site of the diffusion CA) of glucose. For the V79 cell line, a value of 40 104 nmol/cell/h has been used in G1 phase for Vm.
ARTICLE IN PRESS M. Richard et al. / Applied Radiation and Isotopes 67 (2009) 443–446
445
Table 2 Parameters of the survival curve model
Glioma G1 S G2/M V79
as
ar
b
Mean
Standard deviation
0.615 6.803 0.329 0.7
0.03 0.266 0.048 0.1
0.038 0.012 0.043 0.012
0.508 0.113 0.429 0.27
0.33 0.264 0.03 0.12
Mean and standard deviation refer to the threshold dose triggering repair processes (dc, see main text).
characteristic of an individual and was normally distributed within the population of cell. The survival probability was calculated by SFðdÞ ¼ expððas þ b dÞ dÞ
if dodc
SFðdÞ ¼ expððar þ b dÞ dÞ
if dXdc
(6)
In Eq. (6) as4ar are the low-dose constants of the conventional linear-quadratic model, and b is the high-dose constant. The parameters (as, ar, b, mean and standard deviation of dc) should vary throughout the cell cycle. However, in the absence of available data for the simulations of the V79 cell line, the same parameters were used for all the cell phases (see Table 2). For any irradiated cell, the value dc is computed, and the survival probability is determined using the equation relevant for the domain the irradiation dose belongs to. The probability that the hit cell would release a signal was dependent on dose (Schettino et al., 2005): pBS ¼ 4 dose
Fig. 3. Glioma colony at day 72. A single glioma cell was seeded in the middle of the culture dish in a uniform field of nutrient, and allowed to grow. The nutrient was not replenished. At day 72, the colony shows the well-known structure of spheroids: a necrotic core, surrounded by non-proliferative cells and an outermost rim of proliferative cells (dashed lines have been added to help distinguishing the three regions). The radius of the colony is 660 mm and the radius of the necrotic core is about 150 mm.
(7)
If the bystander signal was released, bystander cells had a probability of 88% to survive it (Schettino et al., 2005). The bystander signal could affect any cell on the dish (Schettino et al., 2003, 2005; Prise et al., 1998).
SF
3. Results 3.1. Validation of the model: un-irradiated case Some simulations were first run to test the capabilities of the program to model the growth of a colony in an initially uniform field of glucose. The fate of single glioma cell seeded in the middle of a 144 mm2 dish area was computed and the resulting colony after 72 days is displayed in Fig. 3. The colony radius was 660 mm; the characteristic structure of a spheroid was well reproduced by the model: a necrotic core of dead cells, surrounded by non-cycling cells, and an outermost rim of proliferative cells (Hajikarim and Carlsson, 1978). The necrotic core appeared when the colony had a diameter of about 570 mm, which is consistent with the experimental results from Hajikarim and Carlsson (1978). 3.2. Targeted irradiation of Chinese hamster cells On a 36 mm2 dish, 160 cells were randomly seeded. At the beginning of the simulation, one or all the cells of the dish were given a dose of 0, 0.1 or 0.3 mGy of X-ray irradiation. The simulation was run for a time equivalent to 3 days of culture, and the number of clusters of cells was counted. The survival fraction was calculated as number of cell clusters SF ¼ 160
Survival Probability of V79 cells
(8)
1 0.95 0.9 0.85 0.8 0.75 0.7 0.65 0.6 0
0.5
1.0
1.5
2.0
2.5
3.0
dose (Gy) Fig. 4. Survival probability of V79 cells. The survival probability has been generated by Eq. (6) (see main text). Parameters are displayed in Table 2.
Fig. 4 shows the survival probability curve of the population of V79 cells generated by Eq. (6). Fig. 5 shows the mean survival fraction of the different experiments described above. Each value is the mean of 10 independent simulations run with different seedings of the random number generator. The survival fractions obtained when either one or all cells are irradiated are very similar; the BE is significant for doses above 1 Gy, but not at 0.05 Gy. The standard deviation in case only one cell is irradiated is higher than when all cells are irradiated: this may illustrate the fact that the BE is an ‘‘all or nothing’’ process.
4. Discussion The model described in the present paper was capable of accurately computing the growth of colonies, as shown in Fig. 3. The compartmented structure (necrotic core, resting cells and
ARTICLE IN PRESS 446
M. Richard et al. / Applied Radiation and Isotopes 67 (2009) 443–446
killed (i.e. 0.88); it should be tested whether the probability of bystander killing is modulated, for instance by cell cycle-phase or in a mixed population of several cell lines.
Survival fraction after 3 days culture
1
Survival fraction
0.9 0.8
5. Conclusion
0.7 0.6 0.5 0.4 0
0.05
0.1
0.3
0.5
dose (Gy) 1 cell irradiated
A new mathematical model is presented for the effects of targeted radiation on cell growth. The low-dose HRS and the BE are taken into account and some first results show the resulting survival fraction of a population irradiated with low doses. It is believed that the model can help clarifying issues such as the relationship between HRS and BE, and the mechanisms underlying the BEs. Further investigations are needed on the effect of different radiation types, and on the modulation of the BE with cell-cycle distribution.
All cells irradiated
Fig. 5. Survival fraction in irradiated dishes. Left-hand bars: survival fraction in dishes where one cell was irradiated with 0, 0.05, 0.1, 0.3 or 0.5 Gy. Right-hand bars: survival fraction in dishes where all cells were irradiated.
proliferative cells) was well reproduced; we were not able, though, to compute the saturation in colony radius. This may be due to the absence of growth inhibitory factors in our model (Landry et al., 1982); longer simulations are planned to see whether such a limit appears. The validity of the DCA was also checked against the analytical solution in one dimension (data not shown). Fig. 5 shows that the program is capable of taking account of the HRS and BE. In our simulations, the BE could be induced after 0.1 Gy of irradiation, and the difference in survival fraction between dishes where one cell and all cells were irradiated was not significant; therefore, the BE is predominant in this dose domain. Also, the survival fraction at 0.5 Gy is higher than at 0.3 Gy in dishes where all cells were irradiated: this illustrates the induced radio-resistance following HRS. Further investigations are needed for irradiation with doses above 0.5 Gy. To our knowledge, this is the first model of cell growth, including both the HRS and the BE. The existing models of either HRS or BE were not associated with any population model. Here, the behaviours of single cells (cell cycle progression, response to radiation) are computed and the effect on the population is observed. The results concerning the BEs are preliminary results. The parameters used for the consumption of glucose by V79 cells were guessed, since no value was found in the literature. Also, the survival fraction of V79 to radiation should be further investigated, especially its dependency on cell cycle distribution. It was assumed that the probability for an irradiated cell to trigger the BE was proportional to dose (Eq. (7)) (Schettino et al., 2005). It would be interesting to know if this relationship holds for proton irradiation as well, and for any LET. Also, the model does not allow for bystander cells affected by the signal to be secondary sources. This has yet been suggested by Schettino et al. (2003) as the cause of the clustering of affected bystander cells that the authors showed. The clustering of bystander-induced killed cells has not been seen with the simulations, and thus it may be due to the absence of secondary sources. Finally, a bystander cell receiving the signal had a fixed probability to be
Acknowledgment The work is supported by the Marie Curie Research Network project CELLION, Contract number: MRTN-CT-2003-503923. References Bandini, S., Mauri, G., Pavesi, G., Simone, C., 2001. Parallel simulation of reaction–diffusion phenomena in percolation processes—a model based on cellular automata. Future Gener. Comput. Syst. 17, 679–688. Beuling, E.E., van den Heuvel, J.C., Ottengraf, S.P.P., 2000. Diffusion coefficients of metabolites in active biofilms. Biotechnol. Bioeng. 67, 53–60. Brenner, D.J., Little, J.B., Sachs, R.K., 2001. The bystander effect in radiation oncogenesis: II. A quantitative model. Radiat. Res. 155, 402–408. Casciari, J.J., Sotirchos, S.V., Sutherland, R.M., 1988. Glucose diffusivity in multicellular tumor spheroids. Cancer Res. 48, 3905–3909. Faraday, D.B.F., Hayter, P., Kirkby, N.F., 2001. A mathematical model of the cell cycle of a hybridoma cell line. Biochem. Eng. J. 7, 49–68. Hajikarim, M., Carlsson, J., 1978. Proliferation and viability in cellular spheroids of human origin. Cancer Res. 38, 1457–1464. Joiner, M.C., Marples, B., Lambin, P., Short, S.C., Turesson, I., 2001. Low-dose hypersensitivity: current status and possible mechanisms. Int. J. Radiat. Oncol. Biol. Phys. 49, 379–389. Kirkby, N.F., Burnet, N.G., Faraday, D.B.F., 2002. Mathematical modelling of the response of tumour cells to radiotherapy. Nucl. Instrum. Methods Phys. Res. Sect. B Beam Interact. Mater. Atoms 188, 210–215. Landry, J., Freyer, J.P., Sutherland, R.M., 1982. A model for the growth of multicellular spheroids. Cell Tissue Kinet. 15, 585–594. Little, M.P., Filipe, J.A.N., Prise, K.M., Folkard, M., Belyakov, O.V., 2005. A model for radiation-induced bystander effects, with allowance for spatial position and the effects of cell turnover. J. Theor. Biol. 232, 329–338. Morgan, W.F., 2003. Non-targeted and delayed effects of exposure to ionizing radiation: I. radiation-induced genomic instability and bystander effects in vitro. Radiat. Res. 159, 567–580. Nikjoo, H., Khvostunov, I.K., 2003. Biophysical model of the radiation-induced bystander effect. Int. J. Radiat. Biol. 79, 43–52. Prise, K.M., Belyakov, O.V., Folkard, M., Michael, B.D., 1998. Studies of bystander effects in human fibroblasts using a charged particle microbeam. Int. J. Radiat. Biol. 74, 793–798. Prise, K.M., Schettino, G., Folkard, M., Held, K.D., 2005. New insights on cell death from radiation exposure. Lancet Oncol. 6, 520–528. Richard, M., Kirkby, K.J., Webb, R.P., Kirkby, N.F., 2007. Nucl. Instrum. Methods Phys. Res. Sect. B Beam Interact. Mater. Atoms 255, 18–22. Schettino, G., Folkard, M., Prise, K.M., Vojnovic, B., Held, K.D., Michael, B.D., 2003. Low-dose studies of bystander cell killing with targeted soft X-rays. Radiat. Res. 160, 505–511. Schettino, G., Folkard, M., Michael, B.D., Prise, K.M., 2005. Low-dose binary behavior of bystander cell killing after microbeam irradiation of a single cell with focused C-K X-rays. Radiat. Res. 163, 332–336. Wouters, B.G., Skarsgard, L.D., 1994. The response of a human tumor-cell line to low radiation-doses-evidence of enhanced sensitivity. Radiat. Res. 138, S76–S80.
Contents lists available at ScienceDirect
Applied Radiation and Isotopes journal homepage: www.elsevier.com/locate/apradiso
Cellular automaton model of cell response to targeted radiation M. Richard a,b, K.J. Kirkby a, R.P. Webb a, N.F. Kirkby b, a b
Surrey Ion Beam Centre, Advanced Technology Institute, University of Surrey, Guildford, GU2 7XH, UK Centre for Fluids and Systems, School of Engineering, University of Surrey, Guildford, GU2 7XH, UK
abstract It has been shown that the response of cells to low doses of radiation is not linear and cannot be accurately extrapolated from the high dose response. To investigate possible mechanisms involved in the behaviour of cells under very low doses of radiation, a cellular automaton (CA) model was created. The diffusion and consumption of glucose in the culture dish were computed in parallel to the growth of cells. A new model for calculating survival probability was introduced; the communication between targeted and non-targeted cells was also included. Early results on the response of non-confluent cells to targeted irradiation showed the capability of the model to take account for the non-linear response in the low-dose domain. & 2008 Published by Elsevier Ltd.
1. Introduction
2. Model
The accurate measurement of survival of mammalian cells at low doses of radiation in the past 20 years has shown that it could not be extrapolated from the high dose response. Some cell lines show a hyper radio-sensitivity (HRS) to low doses (below 0.5 Gy), followed by an increased radio-resistance between 0.5 and 1 Gy (Joiner et al., 2001); this leads to survival fraction curves that cannot be described by the traditional linear-quadratic model. On the other hand, not only irradiated cells are affected by radiation, but also bystander non-irradiated cells in the neighbourhood; this has been termed as the bystander effect (BE) and was widely investigated using microbeams (Prise et al., 2005; Morgan, 2003). The mechanisms underlying the HRS and the BE are not fully understood yet. Different experimental procedures have been used and different conclusions, sometimes contradictory, have been drawn. Some models have been created to take into account either the HRS (Wouters and Skarsgard, 1994) or the BE (Little et al., 2005; Nikjoo and Khvostunov, 2003; Brenner et al., 2001). However, none of those models attempted to embed both phenomena. The work presented in this paper and in Richard et al. (submitted) includes the HRS and the BE, using the cellular automaton (CA) computing approach; the diffusion of glucose and the cell cycle are fully described. The response of V79 cells to targeted radiation is simulated.
A CA is a lattice of sites whose time evolution depends on its previous states and on interacting rules between a site and its neighbours. In this model, a first CA computed the diffusion of glucose and a second CA the cell population growth. A detailed description of the diffusion procedure and cell cycle model can be found in Richard et al. (submitted). 2.1. The diffusion of glucose The lattice represented the cellular environment (i.e. culture dish) in which the glucose was diffusing. The diffusion procedure, adapted from Bandini et al. (2001), consisted of two steps. Any site of the CA was divided into four (north, south, east and west) and glucose first diffused between two adjacent quarters of two neighbouring sites (see Fig. 1), following the equation xðt nþ1 Þ ¼ xðt n Þ ðxðt n Þ yðt n ÞÞ k;
D¼
Corresponding author. Tel.: +44 1483686577; fax: +44 1483686581.
E-mail address: [email protected] (N.F. Kirkby). 0969-8043/$ - see front matter & 2008 Published by Elsevier Ltd. doi:10.1016/j.apradiso.2008.06.044
kp0:5
(1)
In Eq. (1), x is the concentration of glucose in one area of a given site (e.g. south) and y the concentration of glucose in the adjacent area of the neighbouring site (e.g. north area of the south neighbouring site).The concentration of glucose is then homogenised on a single site of the lattice. If Dx is the space increment of the CA and Dt the time increment, and k the constant in Eq. (1), the diffusivity D of the glucose is given by
Dx2 k Dt
(2)
Growing cells consumed the diffusing glucose, so its concentration on the lattice decreased. The diffusivity of glucose was chosen to be 7.2 1010 m2/s on the surfaces of the dish free from cells
ARTICLE IN PRESS 444
M. Richard et al. / Applied Radiation and Isotopes 67 (2009) 443–446
(Beuling et al., 2000), and 3 1011 m2/s on the surfaces covered by cells (Casciari et al., 1988).
3
4 x(tn)=1
3
y(tn)=3
2
2
2
5 x(tn+1)=2 y(tn+1)=2
3 1
2.2. The cell cycle
4
The sites of the second CA represented biological cells and the CA evolution showed the growth of the population. The cell cycle was divided into mitosis (M phase) and interphase, i.e. first gap (G1 phase), synthesis (S phase) and second gap (G2 phase) (see Fig. 2). The lengths of M, S and G2 were fixed, while G1 had a probabilistic duration. A G1 cell was allowed to enter the S phase only if the amount of glucose it had absorbed since last division was greater than a given minimum (see Kirkby et al., 2002; Faraday et al., 2001). The consumption of glucose was computed by a Monod-like equation:
5 4
r¼
2 2
2
5 4
1
(3)
When the surrounding density of cells became too high or the surrounding concentration of glucose too low, the cells enter a quiescent state called G]. Values of the different parameters of the model are displayed in Table 1.
5
3
V m CðtÞ K m þ CðtÞ
Fig. 1. Diffusion process. Any site of the CA is divided into four regions, north, south, east and west adjacent to the corresponding neighbouring sites. The glucose first diffuses between the regions of a site and the adjacent regions of the neighbouring sites. The equation of diffusion is x(tn+1) ¼ x(tn)(x(tn))(y(tn)) k with k ¼ 0.5. The glucose concentration is averaged between the four regions of each site.
2.3. Time and space scales of the cellular automata The time scale of the glucose diffusion (order of second) was much smaller than the time scale of the evolution of the cell population (order of hour). As a consequence, it was not possible to choose the same time and space increment for both cellular automata. Therefore, the following relationship was decided:
DxDCA ¼ 5 DxBCA
In Eq. (4), DxDCA is the space increment of the diffusion CA (DCA) and DxBCA the space increment of the biological CA (BCA), which is equal to 12 mm. The time increment of the DCA was then deduced from Eq. (2) to be 0.625 s and the time increment of the BCA was taken to be 1 h. At any iteration of the program, the probabilities of updating either the BCA or the DCA were
G# if not enough nutrient or high cell density Radiation induced cell cycle arrest
M phase G1 phase
Radiation induced cell cycle arrest
G2 phase
pBCA Dt DCA ¼ pDCA Dt BCA
Check nutrient absorption before S Death if absorption
S phase
(4)
Fig. 2. The cell cycle. The cell cycle is divided into mitosis (M phase), during which cells divide, and interphase, during which cells prepare themselves for the next division. The interphase consists of gap 1 (G1 phase), synthesis of DNA (S phase) and gap 2 (G2 phase). The M phase is short compared to the interphase. Cells in G1 continuously check the amount of glucose absorbed since last division; they can enter S only if this amount is greater than an internal threshold. If cells fail to check this criterion after a maximum period of time (maximum length of G1), they die. Cells in G1 can also enter a non-cycling phase, G], if the cell density is too high or if the concentration of glucose is too low. Irradiated cells are arrested in G2, for DNA integrity checking.
(5)
In Eq. (5), pBCA and pDCA are the probabilities of updating the BCA and the DCA, respectively, and DtBCA and DtDCA are the time increments of the BCA and DCA, respectively. 2.4. Microbeam irradiation It was assumed that a cell was very sensitive to a low dose of radiation because the repair processes were not triggered; as the dose increased, those processes would be triggered and the cell would become less sensitive. The limit between the low- and high-dose domain, namely, the dose dc, was supposed to be
Table 1 Main parameters used in the two CA
Glioma V79
DxDCA (mm)
DxBCA (mm)
DtDCA (s)
DtBCA (h)
k
Vm (104 nmol/cell/h)
Km (104 nmol/site)
C0 (104 nmol/site)
60 60
12 12
0.625 0.625
1 1
0.5 0.5
4 40 (G1 phase) 4 (other phases)
75 75
400 400
DxDCA and DxBCA (resp. DtDCA and DtBCA) refer to the space (resp. time) increment of the diffusion and biological CA. k is the constant in Eq. (1) in the main text and C0 is the initial concentration (in mol per site of the diffusion CA) of glucose. For the V79 cell line, a value of 40 104 nmol/cell/h has been used in G1 phase for Vm.
ARTICLE IN PRESS M. Richard et al. / Applied Radiation and Isotopes 67 (2009) 443–446
445
Table 2 Parameters of the survival curve model
Glioma G1 S G2/M V79
as
ar
b
Mean
Standard deviation
0.615 6.803 0.329 0.7
0.03 0.266 0.048 0.1
0.038 0.012 0.043 0.012
0.508 0.113 0.429 0.27
0.33 0.264 0.03 0.12
Mean and standard deviation refer to the threshold dose triggering repair processes (dc, see main text).
characteristic of an individual and was normally distributed within the population of cell. The survival probability was calculated by SFðdÞ ¼ expððas þ b dÞ dÞ
if dodc
SFðdÞ ¼ expððar þ b dÞ dÞ
if dXdc
(6)
In Eq. (6) as4ar are the low-dose constants of the conventional linear-quadratic model, and b is the high-dose constant. The parameters (as, ar, b, mean and standard deviation of dc) should vary throughout the cell cycle. However, in the absence of available data for the simulations of the V79 cell line, the same parameters were used for all the cell phases (see Table 2). For any irradiated cell, the value dc is computed, and the survival probability is determined using the equation relevant for the domain the irradiation dose belongs to. The probability that the hit cell would release a signal was dependent on dose (Schettino et al., 2005): pBS ¼ 4 dose
Fig. 3. Glioma colony at day 72. A single glioma cell was seeded in the middle of the culture dish in a uniform field of nutrient, and allowed to grow. The nutrient was not replenished. At day 72, the colony shows the well-known structure of spheroids: a necrotic core, surrounded by non-proliferative cells and an outermost rim of proliferative cells (dashed lines have been added to help distinguishing the three regions). The radius of the colony is 660 mm and the radius of the necrotic core is about 150 mm.
(7)
If the bystander signal was released, bystander cells had a probability of 88% to survive it (Schettino et al., 2005). The bystander signal could affect any cell on the dish (Schettino et al., 2003, 2005; Prise et al., 1998).
SF
3. Results 3.1. Validation of the model: un-irradiated case Some simulations were first run to test the capabilities of the program to model the growth of a colony in an initially uniform field of glucose. The fate of single glioma cell seeded in the middle of a 144 mm2 dish area was computed and the resulting colony after 72 days is displayed in Fig. 3. The colony radius was 660 mm; the characteristic structure of a spheroid was well reproduced by the model: a necrotic core of dead cells, surrounded by non-cycling cells, and an outermost rim of proliferative cells (Hajikarim and Carlsson, 1978). The necrotic core appeared when the colony had a diameter of about 570 mm, which is consistent with the experimental results from Hajikarim and Carlsson (1978). 3.2. Targeted irradiation of Chinese hamster cells On a 36 mm2 dish, 160 cells were randomly seeded. At the beginning of the simulation, one or all the cells of the dish were given a dose of 0, 0.1 or 0.3 mGy of X-ray irradiation. The simulation was run for a time equivalent to 3 days of culture, and the number of clusters of cells was counted. The survival fraction was calculated as number of cell clusters SF ¼ 160
Survival Probability of V79 cells
(8)
1 0.95 0.9 0.85 0.8 0.75 0.7 0.65 0.6 0
0.5
1.0
1.5
2.0
2.5
3.0
dose (Gy) Fig. 4. Survival probability of V79 cells. The survival probability has been generated by Eq. (6) (see main text). Parameters are displayed in Table 2.
Fig. 4 shows the survival probability curve of the population of V79 cells generated by Eq. (6). Fig. 5 shows the mean survival fraction of the different experiments described above. Each value is the mean of 10 independent simulations run with different seedings of the random number generator. The survival fractions obtained when either one or all cells are irradiated are very similar; the BE is significant for doses above 1 Gy, but not at 0.05 Gy. The standard deviation in case only one cell is irradiated is higher than when all cells are irradiated: this may illustrate the fact that the BE is an ‘‘all or nothing’’ process.
4. Discussion The model described in the present paper was capable of accurately computing the growth of colonies, as shown in Fig. 3. The compartmented structure (necrotic core, resting cells and
ARTICLE IN PRESS 446
M. Richard et al. / Applied Radiation and Isotopes 67 (2009) 443–446
killed (i.e. 0.88); it should be tested whether the probability of bystander killing is modulated, for instance by cell cycle-phase or in a mixed population of several cell lines.
Survival fraction after 3 days culture
1
Survival fraction
0.9 0.8
5. Conclusion
0.7 0.6 0.5 0.4 0
0.05
0.1
0.3
0.5
dose (Gy) 1 cell irradiated
A new mathematical model is presented for the effects of targeted radiation on cell growth. The low-dose HRS and the BE are taken into account and some first results show the resulting survival fraction of a population irradiated with low doses. It is believed that the model can help clarifying issues such as the relationship between HRS and BE, and the mechanisms underlying the BEs. Further investigations are needed on the effect of different radiation types, and on the modulation of the BE with cell-cycle distribution.
All cells irradiated
Fig. 5. Survival fraction in irradiated dishes. Left-hand bars: survival fraction in dishes where one cell was irradiated with 0, 0.05, 0.1, 0.3 or 0.5 Gy. Right-hand bars: survival fraction in dishes where all cells were irradiated.
proliferative cells) was well reproduced; we were not able, though, to compute the saturation in colony radius. This may be due to the absence of growth inhibitory factors in our model (Landry et al., 1982); longer simulations are planned to see whether such a limit appears. The validity of the DCA was also checked against the analytical solution in one dimension (data not shown). Fig. 5 shows that the program is capable of taking account of the HRS and BE. In our simulations, the BE could be induced after 0.1 Gy of irradiation, and the difference in survival fraction between dishes where one cell and all cells were irradiated was not significant; therefore, the BE is predominant in this dose domain. Also, the survival fraction at 0.5 Gy is higher than at 0.3 Gy in dishes where all cells were irradiated: this illustrates the induced radio-resistance following HRS. Further investigations are needed for irradiation with doses above 0.5 Gy. To our knowledge, this is the first model of cell growth, including both the HRS and the BE. The existing models of either HRS or BE were not associated with any population model. Here, the behaviours of single cells (cell cycle progression, response to radiation) are computed and the effect on the population is observed. The results concerning the BEs are preliminary results. The parameters used for the consumption of glucose by V79 cells were guessed, since no value was found in the literature. Also, the survival fraction of V79 to radiation should be further investigated, especially its dependency on cell cycle distribution. It was assumed that the probability for an irradiated cell to trigger the BE was proportional to dose (Eq. (7)) (Schettino et al., 2005). It would be interesting to know if this relationship holds for proton irradiation as well, and for any LET. Also, the model does not allow for bystander cells affected by the signal to be secondary sources. This has yet been suggested by Schettino et al. (2003) as the cause of the clustering of affected bystander cells that the authors showed. The clustering of bystander-induced killed cells has not been seen with the simulations, and thus it may be due to the absence of secondary sources. Finally, a bystander cell receiving the signal had a fixed probability to be
Acknowledgment The work is supported by the Marie Curie Research Network project CELLION, Contract number: MRTN-CT-2003-503923. References Bandini, S., Mauri, G., Pavesi, G., Simone, C., 2001. Parallel simulation of reaction–diffusion phenomena in percolation processes—a model based on cellular automata. Future Gener. Comput. Syst. 17, 679–688. Beuling, E.E., van den Heuvel, J.C., Ottengraf, S.P.P., 2000. Diffusion coefficients of metabolites in active biofilms. Biotechnol. Bioeng. 67, 53–60. Brenner, D.J., Little, J.B., Sachs, R.K., 2001. The bystander effect in radiation oncogenesis: II. A quantitative model. Radiat. Res. 155, 402–408. Casciari, J.J., Sotirchos, S.V., Sutherland, R.M., 1988. Glucose diffusivity in multicellular tumor spheroids. Cancer Res. 48, 3905–3909. Faraday, D.B.F., Hayter, P., Kirkby, N.F., 2001. A mathematical model of the cell cycle of a hybridoma cell line. Biochem. Eng. J. 7, 49–68. Hajikarim, M., Carlsson, J., 1978. Proliferation and viability in cellular spheroids of human origin. Cancer Res. 38, 1457–1464. Joiner, M.C., Marples, B., Lambin, P., Short, S.C., Turesson, I., 2001. Low-dose hypersensitivity: current status and possible mechanisms. Int. J. Radiat. Oncol. Biol. Phys. 49, 379–389. Kirkby, N.F., Burnet, N.G., Faraday, D.B.F., 2002. Mathematical modelling of the response of tumour cells to radiotherapy. Nucl. Instrum. Methods Phys. Res. Sect. B Beam Interact. Mater. Atoms 188, 210–215. Landry, J., Freyer, J.P., Sutherland, R.M., 1982. A model for the growth of multicellular spheroids. Cell Tissue Kinet. 15, 585–594. Little, M.P., Filipe, J.A.N., Prise, K.M., Folkard, M., Belyakov, O.V., 2005. A model for radiation-induced bystander effects, with allowance for spatial position and the effects of cell turnover. J. Theor. Biol. 232, 329–338. Morgan, W.F., 2003. Non-targeted and delayed effects of exposure to ionizing radiation: I. radiation-induced genomic instability and bystander effects in vitro. Radiat. Res. 159, 567–580. Nikjoo, H., Khvostunov, I.K., 2003. Biophysical model of the radiation-induced bystander effect. Int. J. Radiat. Biol. 79, 43–52. Prise, K.M., Belyakov, O.V., Folkard, M., Michael, B.D., 1998. Studies of bystander effects in human fibroblasts using a charged particle microbeam. Int. J. Radiat. Biol. 74, 793–798. Prise, K.M., Schettino, G., Folkard, M., Held, K.D., 2005. New insights on cell death from radiation exposure. Lancet Oncol. 6, 520–528. Richard, M., Kirkby, K.J., Webb, R.P., Kirkby, N.F., 2007. Nucl. Instrum. Methods Phys. Res. Sect. B Beam Interact. Mater. Atoms 255, 18–22. Schettino, G., Folkard, M., Prise, K.M., Vojnovic, B., Held, K.D., Michael, B.D., 2003. Low-dose studies of bystander cell killing with targeted soft X-rays. Radiat. Res. 160, 505–511. Schettino, G., Folkard, M., Michael, B.D., Prise, K.M., 2005. Low-dose binary behavior of bystander cell killing after microbeam irradiation of a single cell with focused C-K X-rays. Radiat. Res. 163, 332–336. Wouters, B.G., Skarsgard, L.D., 1994. The response of a human tumor-cell line to low radiation-doses-evidence of enhanced sensitivity. Radiat. Res. 138, S76–S80.