Biophysical Modeling of the Ionizing Radiation Influence on Cells Using the Stochastic (Monte Carlo) and Deterministic (Analytical) Approaches

Abstract

This review article describes our simplified biophysical model for the response of a group of cells to ionizing radiation. The model, which is a product of 10 years of studies, acts as (a) a comprehensive stochastic approach based on the Monte Carlo simulation with a probability tree and (b) the thereof derived detailed deterministic models describing the selected biophysical and radiobiological phenomena in an analytical manner. Specifically, the presented model describes effects such as the risk of neoplastic transformation of cells relative to the absorbed radiation dose, the dynamics of tumor development, the priming dose effect (also called the Raper–Yonezawa effect) based on the introduced adaptive response approach, and the bystander effect. The model is also modifiable depending on users’ potential needs.

Keywords

low-dose modeling radiation risk adaptive response Monte Carlo dose response low dose

Introduction

From a purely physical perspective, every living organism is a complex system. Such systems are characterized by a very unpredictable behavior, as many different influencing factors must be taken into account. Therefore, it is very difficult, if not impossible, to describe the response of a complex system to a given factor in a strict mathematical formalism. However, using appropriate numerical methods, one can create a model which provides a good approximation of reality within an acceptable range. Such a model can successfully describe the behavior of an organism exposed to one specific external stressor. One of the most important effects of stressor interaction is cancer induction, which is of crucial interest to the presented paper.

Cancer is a genetic disease—it is caused by stable changes in DNA which controls the normal functioning of a cell, created as a result of error or damage. In 2020, 19.3 million estimated worldwide cancer cases were reported, of which almost 10.0 million were fatal.¹ Exposure to external agents may potentially result in induction of a neoplasm in the human body. According to the mutational theories, more than one oncogenic mutation is responsible for initial steps in the neoplastic process, until the cells eventually sustain enough genetic alteration to become autonomous, which results in cancer. It is conservatively estimated that 5% of human cancers are caused by viruses, 5% by radiation, and the remaining 90% by chemicals.²

In the context of the presented paper, ionizing radiation is a special case of an external factor (the stressor) that affects the complex system of a living organism. The subject of the general impact of radiation on health is very wide and will not be discussed here—this article focuses on the description of various attempts to model the impact of ionizing radiation on carcinogenesis initiation and progression, methods of model’s programming and statistical description, as well as on the description of selected biophysical and radiobiological phenomena, which contribute significantly to the complexity of its response to radiation.

Numerous models describe the ionizing radiation interaction with living organism(s), especially on the cellular level. The analysis of different numerical models (Table 1) shows that there is a need for a versatile, adjustable, and homogeneous model that reliably describes the consequences of irradiation: for example, neoplastic transformation (cancer formation) and all significant non-targeted effects. Previous models feature basic consideration of the cellular effects (i.e., bystander effect), are complicated numerical programs using many physical models, or explicitly focus on one type of radiation or cell.

Table 1.

Comparison of different computational stochastic models describing ionizing radiation interaction with the organism on the cellular level.

Paper title/the scope of the model	Reference	Description	Modeling technique	Additional information
“Monte Carlo predictions of DNA fragment-size distributions for large sizes after HZE particle irradiation”	³	DNA damage due to space irradiation	DNA break algorithm (Monte Carlo)	Damage model of large, double-strand DNA due to space radiation
“A biological-based model that links genomic instability, bystander effects, and adaptive response”	⁴	Genomic instability, bystander effect, and adaptive response	Biological-based NEOTRANS₃	For low dose, radiation-induced stochastic effects. Unclear if relates to cancer induction
“A review: Development of a microdose model for analysis of adaptive response and bystander dose response behaviour”	⁵	Dose–response behaviour of radio-protective adaptive response effect	Microdose model	A complex but versatile model that can be applied to many radiation studies. Includes both targeted and non-targeted effects
“Cellular automaton model of cell response to targeted radiation”	⁶	Radiation effect on cell colonies	Cellular automation computing approach	Cell growth model showing the effect of targeted radiation, including both hyper radiosensitivity and bystander effect
“Modeling of radiation-induced bystander effect using Monte Carlo methods”	⁷	Bystander effect	Monte Carlo	The model describes the bystander effect only
“A new view of radiation-induced cancer”	^8-10	Radiation-induced cancer	Mixed deterministic–stochastic formalism	Short-term processes are determined stochastically, while long-term processes deterministically
“A computational model of cellular response to modulated radiation fields”	¹¹	Response of cell populations to spatially modulated radiation exposures	Monte Carlo	Radiation damage is modeled stochastically. The model highlights the influence of intercellular communication (bystander effect). However, this part of the model requires development
“Overview of the PHITS code and its application to medical physics”	¹²	Particle transport	PHITS (Monte Carlo)	General purpose 3D particle transport code. Uses several physical models. Can model all particles over a wide energy range. Does not model biological reactions such as the generation of DNA damage
“The FLUKA code: An accurate simulation tool for particle therapy”	¹³	Interaction of hadrons, heavy ions, and electromagnetic particles	FLUKA code (Monte Carlo)	General purpose code simulating the interaction of hadrons, heavy ions, and electromagnetic particles. Fully integrated particle physics simulation package, using many physical models
“BIANCA, a biophysical model of cell survival and chromosome damage by protons, C-ions and He-ions at energies and doses used in hadrontherapy”	¹⁴	Radiation-induced chromosomal aberrations and resulting in cell death	BIANCA II (Monte Carlo)	Limited radiation types and cell lines. Only one cell death form was described. No mention of non-targeted effects
“TOPAS-nBio: A Monte Carlo simulation toolkit for cell-scale radiation effects”	^15,16	Best for dense track structure and low-fluence scenarios	Track-structure–based Monte Carlo simulations	Extension of the TOPAS Monte Carlo toolkit
“Quantitative modeling of carcinogenesis induced by single beams or mixtures of space radiations using targeted and non-targeted effects”	¹⁷	Space radiation carcinogenesis	Sequential quadratic programming algorithm (Maple-2020 software)	Best-fit dose responses for different radiation types. Focuses on the influence of space radiation. Data based on carcinogenesis of mice. The Monte Carlo procedure is used to estimate uncertainties
“A mathematical radiobiological model (MRM) to predict complex DNA damage and cell survival for ionizing particle radiations of varying quality”	¹⁸	Cell survival after ion irradiation	Deterministic multi-scale approach (MSA)	The radiobiologically based model describes complex DNA damages and cell survival for ion irradiation using three input parameters

This review article shows our own approach which presents selected methods of computer modeling of the effects of ionizing radiation on the body at the cellular level. The authors present a stochastic model using the Monte Carlo technique based on a Markov probability tree with memory approach, which can be successfully used for computer simulation of various types of radiation phenomena. A deterministic approach is also described, which focuses on specific processes, such as tumor growth, or the phenomenon of radiation adaptive response.

Monte Carlo Stochastic Modeling of Irradiated Cells

The Monte Carlo modeling technique has been implemented in various types of biophysical simulations for many years. The development of the model discussed in this chapter started 10 years ago¹⁹ when the first simple Monte Carlo Markov Chain (MCMC) with probability tree model of cellular radiation response was published. Over the years, this model has been refined, bringing its results closer to reality.^20–25 Details of the current version of the model are described below.

Cell Arrangement

The model can simulate the behavior of any number of cells (N), limited only by computer capabilities. All cells are organized into a three-dimensional matrix that simulates real tissue (see Figure 1 and 2A and B). This matrix may undergo changes (e.g., cells may die, multiply, and mutate) with successive discrete time steps. Moving from one time step to the next means that each cell has passed through the probability tree (described in detail later). At each time step, each cell has one of the following states: (a) healthy, (b) damaged, (c) mutated, (d) cancerous, and (e) dead cell/lack of cell (empty). Of course, all those states are on the general level of biological processes.

Figure 1.

Cell’s numerical organization in the three-dimensional matrix which geometry is presented in Figure 2.

Figure 2.

Two possibilities of geometrical cell arrangements in the 3D matrix: (a) the cube or (b) sphere, where 3 axes represent geometric coordinates. The plot (c) shows the schematic possibility of radiation beam implementation.²⁵

Damaged cells are those whose DNA has been successfully hit by ionizing radiation: the act of ionization in the space of DNA has occurred (irrespective of the type of lesion or radiation). A mutant cell, on the other hand, is a cell that has transformed from a damaged cell because of damage that was not repaired or was repaired incorrectly (via all potential ways). A cancerous cell is a neoplastically transformed mutant cell, that is, that has undergone the process of cytogenetic transformation due to the accumulation of an appropriate number of oncogenic mutations in its DNA.²⁶ More details about cell statuses and their transitions with time (number of time steps) will be described later.

Dosimetry

In the simplest variant, all cells are pre-assigned with a certain dose of ionizing radiation. This may be a single dose value (D) at a specific fixed time step (k) or a fixed dose value in all subsequent time steps which corresponds to a constant dose-rate. The two scenarios mentioned are the most commonly used, although it is possible to enter an arbitrarily variable dose-rate (dose per step) as well as to differentiate the dose depending on the cells’ matrix region. In particular, by applying microdosimetric methods, it is possible to simulate any dose distribution based on the adopted external radiation field. For example, a dose distribution for tumor cells was simulated for the interaction of alpha radiation (Bragg curve) in a specific direction of the cell matrix.²⁷ Another example is a possibility of cell irradiation by a specific beam with its dedicated dimensions and angle²⁵ (see Figure 1C).

Probability Tree

At each successive time step, each cell is tested against a probability tree (the so-called Bernoulli tree). The iterative algorithm in three loops passes through each i-th cell in turn and first checks its status (healthy/damaged/mutated/cancerous), which determines further treatment. Then, the algorithm uses Monte Carlo methods to check the appropriate probability functions, such as whether the cell was hit by radiation, whether the cell died naturally, whether the cell divided (creating a daughter cell), or if any further biophysical effects occurred. The diagram of the probability tree used in the model is presented in Figures 3 -7. A description of all probability functions is presented in Table 2. More details are presented in References.^20,24,25 However, one has to note that all biophysical processes are implemented on a rather general level, where many details (e.g., chromosomal aberrations) are just a part of single phenomenon (e.g., mutation creation).

Figure 3.

Initial part of the Monte Carlo stochastic tree of probabilities. See Figures 4 -7 for details in each eight scenarios. The exemplary scheme of iteration is presented in Figure 8. The probability tree was based on its previous published version.^19,20

Figure 4.

Part of a probability tree from Figure 3 concerning healthy cell which was (a) hit and (b) not hit by ionizing radiation.

Figure 5.

Part of a probability tree from Figure 3 concerning damaged cell which was (a) hit and (b) not hit by ionizing radiation.

Figure 6.

Part of a probability tree from Figure 3 concerning mutated cell which was (a) hit and (b) not hit by ionizing radiation.

Figure 7.

Part of a probability tree from Figure 3 concerning cancerous cell which was (a) hit and (b) not hit by ionizing radiation.

Table 2.

List of probability functions used in the Monte Carlo model (Figures 3 –7), with their descriptions and parameters.

Equation	Description
$P_{h i t} = 1 - e^{- a D}$	The probability that a cell is hit by radiation, where D is the dose per time step (for a type of radiation) and a is an input constant (analogically in next items)
$P_{D, D D, M D} = (1 - τ) (1 - e^{- a K^{n}}) + τ$	The probability of a natural death of a healthy (D), damaged (DD), or mutated (MD) cell, where K, n, and τ are the cell’s age, shape constant, and scaling constant, respectively
$P_{C D} = c o n s t$	The probability of natural death of a cancerous cell
$P_{R D} = 1 - e^{- a D}$	The probability that a cell dies due to irradiation
$P_{S, D S, M S, C S} = c o n s t$	The probability that a healthy (S), damaged (DS), mutated (MS), or cancerous (CS) cell multiplies; please note, that the new cell (daughter cell) can be located within some r=4 radius around the mother cell
$P_{M} = (1 - τ) (1 - e^{- a K^{n}}) + τ$	The probability of spontaneous damage occurrence in a cell; for simplicity, this relation can be approximated by a linear function, P_M = τ + a K
$P_{R D E M} = 1 - e^{- a D}$	The probability that a cell gets damaged due to irradiation
$P_{R} = q e^{- a K^{n}}$	The probability of natural repair of one of the lesions in a cell (each type), where q corresponds to the experimental constant
$P_{D M} = 1 - e^{- a U}$	The probability of a cell mutation, where U is the number of damages
$P_{R M M} = 1 - e^{- a U_{2}}$	The probability of a new (next) mutation in a mutated cell where $U_{2}$ is the number of lesions accumulated after the previous mutation
$P_{R C} = 1 - e^{- a M^{n}}$	The probability that a mutated cell transforms into a cancerous cell, where M is the number of mutations, a is an empirical constant, and n is a critical index²⁸
$P_{C R D} = 1 - e^{- a D}$	The probability that a cancerous cell dies because of its radiosensitivity
$P_{A} = \sum_{k = 0}^{K} α_{0} D_{k}^{2} {(K - k)}^{2} e^{- α_{1} D_{k} - α_{2} (K - k)}$	The cumulative probability of an adaptive response, where k, D_k, and K are past time step, dose in k step, and cell’s age, respectively; parameters {α} represent experimental constants
$P_{B} (t + Δ t) = P_{B} (t) + λ_{M D P} M Δ t + λ_{G J P} G Δ t$ or $P_{B} = β_{1} (1 - e^{- β_{2} D})$ and $P_{B}^{'} (r) = \frac{β_{3}}{r!}$	The probability of a bystander effect in two different approaches: First one based on the Japanese model^29,30 (see equation (4) and its description), and the second one is a two-step approach, where r represents the distance between the bystander and irradiated cell and {β} are empirical constants
	See equations no. 32 and 33 in Reference 30 for the first approach

It should be emphasized that the passage of all cells in the matrix through the probability tree (each cell individually) equals one time step. In the following time step, the situation repeats itself and again all cells go through the probability tree according to the same pattern (see Figure 8). However, due to the random nature of probability, the results will be different each time. This approach allows for an evolutionary change of the entire cell layout depending on time, dose, and other model input parameters.

Figure 8.

Simplified exemplary scheme of the algorithm application to a cell colony. The colony is composed in the first step of healthy cells only (top line), which moves through the probability tree (Figures 3 -8) to create the next step colony.

Probability Functions

The probability tree described in the previous section is based on a number of probability functions that describe biophysical processes in a cell. Some of them are simple functions, described by constant values, but most of them are complex functions of many variables, which include, among others, time, age of the i-th cell, dose in a given time step, radiation doses in previous time steps, number of lesions, or number of mutations. These variables are accompanied by appropriate constants, the so-called model parameters, the determination of which is non-trivial. These parameters can be different for various types of cells, different tissues, or organisms. Therefore, one of the basic activities is to calibrate the model based on the existing experimental data, allowing determination of the values of parameter inputs to the model for specific conditions. Such calibration is difficult and of high complexity but is also necessary for the correct quantitative operation of the model. Of course, arbitrary values of the input parameters could be assumed; however, the results of such simulations would be only qualitative and would not relate to any particular cell type or experimental case.

As mentioned earlier, the model contains many probability functions that comprehensively, albeit on a general level, describe the biophysics of irradiated cells. A detailed description of these functions, along with their mathematical formulae, is given in Table 2. Some of them describe quite obvious biological processes, such as death, multiplication, or radiation damage repair, but some of them address less known phenomena, called non-targeted effects, which include the radiation adaptive response and the bystander effect. Some of them are described below in a more precise way.

Radiation Adaptive Response

The phenomenon of radiation adaptive response includes stimulating the body with low doses of radiation in order to stimulate its repair mechanisms.^31-35 This phenomenon is very subtle and does not always occur;³³ hence, there is a strong necessity to describe it using stochastic methods. If it does occur, it plays an essential role in the cell response to ionizing radiation in the low-dose area and is one of the factors that generate the so-called non-linear response.^21,36

In the presented model, an equation derived from other studies is used, describing the probability of an adaptive response for the dose D received k time steps ago:^20,21,37

p_{A} = α_{0} D^{2} k^{2} e^{- α_{1} D - α_{2} k}

(1)

where {α} are input parameters (described later). The above equation de facto describes the probability of repairing any damage to the cell’s DNA, not necessarily related to radiation, some time (k) after the dose D. The exemplary time-dependent shape of equation (1) is presented in Figure 9 (after constant single D). As a cell may be irradiated with different doses at different time steps, in practice the sum of the independent equations (1) is used for each dose separately²⁰ (see Table 2).

Figure 9.

Exemplary hunchbacked shape of equation (1) representing the radiation adaptive response probability function, after single constant dose irradiation.

A special case of the radiation adaptive response is the so-called priming dose effect, also known as the Raper–Yonezawa effect.³⁸ It describes the response to a situation where a given organism is first irradiated with a low priming dose and after some time with a high challenging dose. Due to this effect, the negative health effects on the body are smaller than if it was exposed to only a high dose. More about the Raper–Yonezawa effect will be provided later.

Bystander Effect

The bystander effect is a classic example of the non-targeted effect of the interaction of ionizing radiation. It describes the phenomenon where the non-irradiated cell adjacent to the one irradiated behaves as if it were irradiated itself. This is due to the transmission (via diffusion) of molecular signals from a cell hit by radiation to a neighboring cell. Firstly, the probability of the bystander effect occurring depends on the dose-related function, which is saturated for medium and high doses:^20,30,39

P_{B} = β_{1} (1 - e^{- β_{2} D})

(2)

and secondly, after fulfilling the probability given by equation (2), the empirical distribution is given by the Poisson distribution:²¹

P_{B}^{'} (r) = \frac{β_{3}}{r!}

(3)

where r represents the distance between bystander and irradiated cell and {β} are empirical constants. As a result, the bystander cell receives a pulse analogous to a hit by radiation, which is assumed to be detrimental.³⁰

The presented description is not the only approach to the bystander effect. A good alternative is based on the diffusion equation with 2 ways of the signal transmission, analogically as in Japanese model,²⁹ which can be directly implemented in the presented model³⁰ using the simplified equation:

P_{B} (t + Δ t) = P_{B} (t) + λ_{M D P} M Δ t + λ_{G J P} G Δ t

(4)

where Δt is the time interval of bystander signal dispersion (usually — a time step) and λ_MDP and λ_GJP are experimental coefficients related to two intercellular signaling pathways — the medium-mediated pathway (MDP) and the gap junctional pathway (GJP), respectively. The most important parameters are dose- and time-related factors M and G (which are precisely described in equations no. 32 and 33 in Reference 30). Because of their length and additional parameters, the details of those functions will be omitted here.

To conclude, equation (4) was finally used in the presented model with parameters from Table 3.

Table 3.

Estimated values of input parameters for probability functions (see Table 2) of Monte Carlo stochastic model, dedicated for human lymphocytes and their response to ionizing radiation (X-rays). The parameter estimation was obtained by partial model calibration on experimental data for dedicated processes.^3,29,39–51

Effect from Table 2	Value of Parameter
Cell is hit by radiation	$a_{h i t} = 1.3 {G y}^{- 1}$
Natural death of a healthy cell	$τ_{D} = 0.00035$ $a_{D} = 2 \cdot 10^{- 12} h^{- 3}$
Natural death of a healthy cell	$n_{D} = 3$
Natural death of a damaged cell	$τ_{D D} = 0.00045$ $a_{D D} = 5 \cdot 10^{- 12} h^{- 3}$
Natural death of a damaged cell	$n_{D D} = 3$
Natural death of a mutated cell	$τ_{M D} = 0.001$ $a_{M D} = 1 \cdot 10^{- 10} h^{- 3}$
Natural death of a mutated cell	$n_{M D} = 3$
Natural death of a cancerous cell	$P_{C D} = 0.0004$
Death due to irradiation	$a_{R D} = 0.5 G y^{- 1}$
Multiplication of a healthy cell	$P_{S} = 0.0027$
Multiplication of a damaged cell	$P_{D S} = 0.002$
Multiplication of a mutated cell	$P_{M S} = 0.002$
Multiplication of a cancerous cell	$P_{C S} = 0.009$
Spontaneous damage	$τ_{M} = 0.001$ $a_{M} = 6.8 \cdot 10^{- 12} h^{- 3}$
Spontaneous damage	$n_{M} = 3$
Damage due to irradiation	$a_{R D E M} = 2.4 G y^{- 1}$
Natural repair of one lesion	$q_{R} = 0.04$
	$a_{R} = 2.5 \cdot 10^{- 12} h^{- 4}$
	$n_{R} = 4$
Mutation creation in a damaged cell	$a_{D M} = 1 \cdot 10^{- 6}$
Subsequent mutation in a mutated or cancerous cell	$a_{R M M} = 2 \cdot 10^{- 4}$
Neoplastic transformation	$a_{R C} = 2 \cdot 10^{- 9}$
Neoplastic transformation	$n_{R C} = 9$
Death of a cancerous cell due to radiosensitivity	$a_{C R D} = 0.32 G y^{- 1}$
Adaptive response	$α_{0} = 22.9 {G y}^{- 2} h^{- 3}$
	$α_{1} = 79.4 {G y}^{- 1}$
	$α_{2} = 0.0832 h^{- 1}$
Bystander effect—Japanese approach (in Reference 30)	$d = 1 \cdot 10^{- 5} m$
	${}^{M}W = 1 \cdot 10^{- 10} m^{2} h^{- 1}$
	${}^{M}α = 1 {G y}^{- 1}$
	${}^{M}β = 4.6 \cdot 10^{- 6} h^{- 1}$
	${}^{G}{W = 5 \cdot 10^{- 11}} m^{2} h^{- 1}$
	${}^{G}{α = 1} {G y}^{- 1}$
	${}^{G}{β = 1.18 \cdot 10^{- 3} h^{- 1}}$
	$λ_{M D P} = 0.01$
	$λ_{G J P} = 0.01$
	$Δ t = 0.1 h$
Bystander effect—own approach^a	$β_{1} = 0.08$
	$β_{2} = 50 {G y}^{- 1}$
	$β_{3} = 0.2$

^aParameters {β} for bystander effects were calculated for human-hamster hybrid fibroblasts irradiated by α radiation³⁹ (see equations (2) and (3)).

Neoplastic transformation

The neoplastic transformation (cancer transformation) is a key element of the presented model, as it was its original purpose. Due to this, it is possible to evaluate the so-called radiation risk, that is, the overall probability of cancer formation for a specific radiation dose.³⁷

The described model uses the sigmoid function (the same as for the phase transition dynamics) to describe the probability of neoplastic transformation depending on the number of M mutations with critical index n:

P_{R C} = 1 - e^{- a M^{n}}

(5)

The above formula is known as the Avrami equation and shows that only a certain accumulated number of M mutations determines the highly non-linear probability that a mutated cell changes its status to a cancerous one. The derivation of this equation is based on the theory of nucleation and growth of crystals, as an analogue to cancer.²⁸

Applicability

The presented model has wide applications in simulating various cellular processes, but its main goal is to simulate the processes of radiation carcinogenesis. This means that the model can be used to support radiation risk analyses and can be used, for example, for strictly radiobiological purposes, as well as auxiliary in radiological protection or probabilistic safety analyses (PSA) of the Level 3 in a nuclear power plant and its emergency planning.

Within the presented model, any initial conditions can be set (e.g., all damaged or healthy cells, a specific number of cells and their spatial distribution, and specific dose fractionation). Therefore, it is possible to simulate, for example, radiotherapeutic processes, where cells with a set status of “cancer” are exposed to high doses of radiation, and neighboring healthy cells may be affected by, for example, the bystander effect. This is one of many examples of how the model can be used.

In this context, an interesting example of the results that can be achieved with the model is the simulation of the cell survival curve for a group of cells with “cancer” status. Based on data from the in vitro experiment for the U251 human glioblastoma and D384 medulloblastoma cell lines,⁵² the model allows for a quantitative assessment of the effect of cancer cell killing.²⁵ Additionally, the application of the ionizing radiation beam (Figure 2C) enables the simulation of a radiotherapy treatment process by directing the beam to a specified location in the simulation space with an appropriately entered angle as an additional input parameter of the simulation.

What’s more, additional work is ongoing to develop the model’s capabilities and thus increase its applicability for potential users.

Exemplary Calibration

Calibration of the model requires determining the values of parameters used in the probability functions, so that the simulation results match experimental data for exact cell type. In particular, this chapter focuses on the exemplary calibration process on human lymphocytes as the most frequently used type of cells. Parameters from each function were determined individually, without including other mechanisms in the model, unless otherwise stated. Subsequent cellular effects were then included in the simulation to determine values more accurately.²⁴

There were different ways in which the calibration was done depending on the function type. In most cases, the calculated parameter was a mean value that characterizes the entire population, for example, mutation rate, division rate, or a number of lesions per cell. In some cases, two functions were related. For example, the number of damages in a cell depends on the probability that spontaneous damage occurs in a cell and on the probability of natural repair of one of the lesions in a cell. In those cases, parameter values were initially estimated for each individual function and then changed with each subsequent simulation. This allowed the observation of how parameter alterations influenced the obtained results and to identify which changes obtained results similar to those presented in experimental studies.⁵⁰

To determine the parameters describing the probability of cell multiplication and cell death, the mean frequency of cell division in a stable colony (where all the elements of the array are occupied by cells) was used, which was determined as an average for different types of lymphocytes.⁴⁵

The probability of cell death from radiation depends on whether it was hit or not. That is why the parameter used for calculation of the probability of death directly from radiation (a_RD) was determined simultaneously with the parameter used for finding the probability of a cell being hit (a_hit). Survival curves were generated for subsequent values of the parameter a_hit and a_RD, and a chi-square test was performed to see which curve was the best fit to the curve presented in the experimental study.⁴⁴

Parameters associated with the occurrence of mutations in a cell were determined for non-irradiated cells.⁴³ The constants for determining the probability of radiation damage were then adjusted to generate a correspondingly higher number of mutations after irradiation.⁴¹ Then, because the frequency of occurrence of mutated cells is low, parameters for determining the probability of subsequent mutation⁴³ and malignant transformation⁴⁸ were determined using only the mutated cell population.

In the determination of parameters describing tumor growth,⁴⁹ a higher density of cancerous tumors was included by reducing the radius (r) around the mother cell (where a new cell can be located) by a value of 1. The value of the parameter used to determine the probability of death of a cancer cell due to its radiosensitivity was determined by finding the best fit of the survival curve—obtained in the model—to the experimental curve⁴⁰ while considering cancer cells that died due to irradiation.

Due to limitations in the available experimental data, calibration of {β} parameters (equations (2) and (3)) related to the bystander effect was carried out for fibroblasts irradiated with alpha particles. The number of damages induced by the bystander effect dominates for doses <0.2 Gy; therefore, the maximum value of the bystander effect probability was determined as the intersection point with the function describing the probability of radiation damage for a dose equal to 0.2 Gy. The parameter determining the reach of the bystander effect was chosen so that for a colony of 100 cells approximately 10% of them were damaged.³⁹

Parallelly, the values of the parameters used to describe the bystander effect according to the Japanese model (equation (4)) were taken directly from a paper describing the original model²⁹ (see Table 2). Analogically, parameters {α} describing the adaptive response (see equation (1) and Table 2) were taken directly from Reference 38.

Estimated values of all parameters (related to human lymphocytes) are presented in Table 3. These are calculated values based on the experimental data from specific papers and, due to large discrepancies in available experimental data, do not guarantee the consistency of simulation results with data from other articles, especially for different types of cells.

Deterministic Modeling of Irradiated Cells

Monte Carlo modeling uses stochastic formalism as it is the only possible mathematical description of a complex process such as the response to the irradiation of group of cells. It is not possible to present an equivalent single analytical description of this entire process, only of a selected process; this can be extracted from the probability tree and described with one analytical deterministic mathematical formula. Selected examples of such deterministic approaches, which can be considered as separate models describing specific biophysical phenomena, are presented and discussed below.

“Lesion to Cancer” Model

Based on the previously presented Monte Carlo model, it is possible to analytically present a narrow path from a broad probability tree, starting from the factor inducing a single damage to the neoplastic transformation.³⁷ This approach is called the “Lesion to Cancer” model⁵³—it covers all processes from initial lesion (damage), through repair (or lack thereof), the generation and accumulation of mutations and, consequently, neoplastic transformation. This model focuses on the biophysical impact of a DNA region that is being hit. Based on general knowledge, it was assumed that in order to initiate the carcinogenic process, it is crucial to cause a mutation in the DNA coding region, preferably either in the oncogene or tumor suppressor gene. This is, however, a simplification of this complex mechanism, but it allowed us to make the model more biologically accurate. Thanks to the applied approaches, it was possible to develop one mathematical formula describing the number of cells that underwent neoplastic transformation:³⁷

N (t, M) = N_{0} {(1 - e^{- B P_{M} t})}^{M} (1 - e^{- c M^{k}})

(6)

where t denotes the time after initial irradiation, M is the number of mutations created, N₀ is the initial number of healthy cells, B is the number of critical DNA bases in critical codons of all tumor-associated genes per cell, c is an empirical constant approximately equal to 0.028, k is a critical index for 1 dimensional growth (k = 2), and P_M is a probability of mutation given by formula no. (13) presented in the paper by Dobrzyński et al.³⁷ One can assume, for simplicity, that BP_M ≈ 0.01 year⁻¹. Note that the last term in equation (6) corresponds to the Avrami formula, see equation (5).

Tumor Growth Model

For longer periods of time (t→∞), equation (6) presented in the previous subchapter can be summed over all values of M to obtain the tumor growth function.³⁷ This has been discussed in a recent paper by Fornalski et al,⁵⁴ where the Gompertz function for the dynamics of tumor growth was precisely described based on biophysical grounds. Additionally, the herein presented calculations have been independently supported by the Monte Carlo simulations^23,25 demonstrating that the Gompertz curve does describe cancer cells dynamics adequately, except at the very early phase where a parabolic function is more appropriate. These results are consistent with the previous reports by other authors.^55,56

Tumor growth modeling is a topic which can be discussed in the context of the dose, time, or even dose-rate relationships. In this regard, the Gompertz curve is not the only model that can be considered: other sigmoidal functions can, to a certain extent, also properly describe the tumor dynamics.³⁶

The Priming Dose Effect Model

The recent example of the deterministic solution related to the adaptive response description only is the one dedicated to the priming dose effect (also called the Raper–Yonezawa effect, which was mentioned earlier). In this case, the radiation adaptive response was taken as a theoretical background to calculate the percentage decrease (δ) of post-radiation mutations after the irradiation scheme of priming + challenging doses (D₁+D₂) in relation to the isolated challenging dose (D₂) case:³⁸

δ = 1 - e^{- ξ_{D 1} e^{- α_{2} ∆ t} [α_{2} ∆ t (α_{2} ∆ t + 2) + 2]} - \frac{D_{1}}{D_{2}} e^{- 2 ξ_{D 1}}

(7)

where

ξ_{D_{1}} = \frac{α_{0}}{α_{2}^{3}} D_{1}^{2} e^{- α_{1} D_{1}}

(8)

and {α} are mentioned earlier model’s free parameters (obtained via its calibration), while Δt corresponds to the time interval between small priming dose (D₁) and the large challenging dose (D₂). The diagram of the Raper–Yonezawa effect is presented in Figure 10.

Figure 10.

Diagram of the Raper–Yonezawa effect (also known as the priming dose effect): the single priming (small) dose (D₁) generates much fewer negative effects than the single detrimental (high) dose (D₂). However, when D₂ follows D₁ with some time distance between them (Δt), the overall biological effect for that D₁+D₂ total dose is lower than for single D₂. The parameter δ is therefore showing the percentage difference between the effect (e.g., mutations or lesions) generated by the single dose D₂ (without the priming dose) and the combination of D₁+D₂.

The presented model works well to describe the practical application of the phenomenon of radiation adaptive response. It was calibrated on various experimental data, such as lesions in human lymphocytes and chromosomal inversions in mice.³⁸

Constant High Background Radiation Model

Another popular aspect of the implementation of radiation adaptive response is the case of high background radiation areas (HBRAs), where some epidemiological studies show significant radioadaptation.^57,58 However, the presented theoretical approach given by equation (1) is related to the single-dose pulse (D) while HBRA corresponds to the constant dose-rate ( $\dot{D}$ ) irradiation with continuous (not discrete) time (k→t). The only way to modify this approach is to sum all individual adaptive response signals given by equation (1) which results in sigmoidal (namely, saturated for long time) probability function of⁵⁹

P_{C} = \int_{t = 0}^{\infty} P_{A} (\dot{D}, t) d t = μ_{0} {\dot{D}}^{2} e^{- μ_{1} \dot{D}}

(9)

where {µ} are empirical constants analogical to {α} but related to the dose-rate. For practical purposes, it is easier to use the repair effectiveness, R = 1 – exp(-P_C),⁵⁹ which can be directly related to epidemiological studies of individuals exposed to elevated values of constant dose-rates.

Discussion and Conclusions

This review article describes our proposed stochastic model of a cell group response to ionizing radiation using the Monte Carlo technique with a probability tree. This universal but simplified approach allows for modification and the creation of virtually any model: depending on the needs, the appropriate branches of the tree can be expanded; thus, enabling a more detailed simulation of the selected phenomenon. This model can be calibrated (i.e., its input parameters estimated) based on existing or dedicated experimental data describing a specific cell type, both human and animal.

Thanks to the applied stochastic approach, it is possible to treat a set of cells as a complex physical system, the precise analytical description of which seems to be currently impossible. However, using computer simulations, basically any simplified behavior of such cells can be modeled. An interesting case in particular is the determination of the probability of neoplastic transformation for a given dose, which can be used to determine the radiation risk (the risk of stochastic effects), especially cancer risk. For such an estimate to be reliable, the model takes into account many different types of biophysical phenomena, such as repair mechanisms, as well as the non-targeted effects, namely, the bystander effect and the adaptive response.

The presented model can be simplified and narrowed in certain ranges and, for example, only a selected phenomenon can be described, which can be demonstrated in the form of one narrow branch in a broad probability tree. Such a narrow branch, with the probability functions included within it, can therefore be written in an analytical form, creating de facto a new (secondary) deterministic model, which is, in fact, an independent fragment of the main model. So far, five such secondary deterministic models have been developed,^{36–38,54,59} but more are in preparation. This means that the presented approach is broad and universal which enables a relatively accurate biophysical description of various aspects of irradiated cells within natural model limits. Although each model is just an approximation of the reality, its proper implementation with validated input data allows for irradiated cellular studies, including cancer treatment, prognosis, or even prevention.

Footnotes

Acknowledgments

The authors would like to kindly thank Prof. Marek Janiak (Poland) and Dr Yehoshua Socol (Israel) for substantial help during the model development. Additional thanks to Prof. Marcin Kruszewski (Poland) and Dr Sylwester Sommer (Poland) for consultations in radiobiology and to Ms Charley Jeynes (University of Manchester) for linguistic corrections.

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The research was funded by the National Centre for Nuclear Research (NCBJ), Poland.

ORCID iDs

Krzysztof W. Fornalski

Łukasz Adamowski

References

Sung

Ferlay

Siegel

Laversanne

Soerjomataram

Jemal

Bray

. Global Cancer Statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 Cancers in 185 Countries. CA A Cancer J Clin. 2021;71(3): 209-249. doi:10.3322/caac.21660.

Malarkey

Hoenerhoff

Maronpot

. Carcinogenesis: Mechanisms and Manifestations. Haschek and Rousseaux's Handbook of Toxicologic Pathology (Third Edition). 2013;I, Chapter 5:107-146. doi:10.1016/B978-0-12-415759-0.00005-4.

Ponomarev

Cucinotta

Sachs

Brenner

. Monte Carlo Predictions of DNA Fragment-Size Distributions for Large Sizes after HZE Particle Irradiation. Phys Med. 2001;17(suppl 1). In:Proceedings of the 1st International Workshop on Space Radiation Research and 11th Annual NASA Space Radiation Health Investigators’ Workshop, Arona (Italy), May 27-31, 2000.

Scott

. A biological-based model that links genomic instability, bystander effects, and adaptive response. Mutat Res Fund Mol Mech Mutagen. 2004;568(1):129-143.

Leonard

. A review: Development of a microdose model for analysis of adaptive response and bystander dose response behavior. Dose-Response. 2008;6:113-183.

Richard

Kirkby

Webb

Kirkby

. Cellular automaton model of cell response to targeted radiation. Appl Radiat Isot. 2009;67(3): 443-446.

Xia

Liu

Xue

Wang

. Modeling of radiation-induced bystander effect using Monte Carlo methods. Nucl Instrum Methods Phys Res Sect B Beam Interact Mater Atoms. 2009;267(6):1015-1018.

Shuryak

Hahnfeldt

Hlatky

Sachs

Brenner

. A new view of radiation-induced cancer: integrating short- and long-term processes. Part I: Approach. Radiat Environ Biophys. 2009;48: 263-274.

Shuryak

Hahnfeldt

Hlatky

Sachs

Brenner

. A new view of radiation-induced cancer: integrating short- and long-term processes. Part II: second cancer risk estimation. Radiat Environ Biophys. 2009;48:275-286.

10.

Shuryak

Sachs

Brenner

. A new view of radiation-induced cancer. Radiat Protect Dosim. 2011;143(2-4):358-364.

11.

McMahon

Butterworth

McGarry

, et al. A Computational model of cellular response to modulated radiation fields. International Journal of Oncology, Biology, Physics. 2012:84(1): 250-256.

12.

Sato

Niita

Matsuda

Hashimoto

, et al. Overview of the PHITS code and its application to medical physics. Progress in Nuclear Science and Technology. 2014;4:879-882.

13.

Battistoni

Bauer

Boehlen

, et al. The FLUKA code: An accurate simulation tool for particle therapy. Front Oncol. 2016. doi:10.3389/fonc.2016.00116

14.

Carante

Aime

Cajiao

JJT

Ballarini

. BIANCA, a biophysical model of cell survival and chromosome damage by protons, C-ions and He-ions at energies and doses used in hadrontherapy. Phys Med Biol. 2018;63(7).

15.

Schuemann

McNamara

Ramos

, et al. TOPAS-nBio: A Monte Carlo simulation toolkit for cell-scale radiation effects. NASA Report; 2019, https://three.jsc.nasa.gov/articles/TOPAS-nBio_Schuemann.pdf

16.

Schuemann

McNamara

Ramos-Mendez

, et al. TOPAS-nBio: An Extension to the TOPAS Simulation Toolkit for Cellular and Sub-cellular Radiobiology. Radiat Res. 2019;191(2):125-138.

17.

Shuryak

Sachs

Brenner

. Quantitative modeling of carcinogenesis induced by single beams or mixtures of space radiations using targeted and non-targeted effects. Sci Rep. 2021;11(1):1-12.

18.

Kalospyros

Nikitaki

Kyriakou

Kokkoris

Emfietzoglou

Georgakilas

. A mathematical radiobiological model (Mrm) to predict complex dna damage and cell survival for ionizing particle radiations of varying quality. Molecules. 2021:26:840.

19.

Fornalski

Dobrzyński

Janiak

. A Stochastic Markov Model of Cellular Response to Radiation. Dose-Response. 2011;9(4):477-496.

20.

Fornalski

. Mechanistic model of the cells irradiation using the stochastic biophysical input. Int J Low Radiat. 2014;9(5/6): 370-395.

21.

Fornalski

Dobrzyński

Reszczyńska

2017. Modelling of the radiation carcinogenesis: the analytic and stochastic approaches. In: Ainsbury

Calle

Cardis

Einbeck

Gómez

Puig

, (eds). Extended Abstracts Fall 2015. Trends in Mathematics, vol. 7, subseries: Research Perspectives CRM Barcelona (Springer); 2017:95-101. doi:10.1007/978-3-319-55639-0_16.

22.

Wysocki

. Modelowanie odpowiedzi grupy komórek na promieniowanie jonizujące metodą Monte Carlo. Warszawa, Poland: Engineer Thesis, Faculty of Physics, Warsaw University of Technology; 2017. (supervisors: K.W. Fornalski, I. Słonecka).

23.

Wysocki

. Modelowanie Monte Carlo efektu sąsiedztwa oraz dynamiki procesu nowotworzenia dla grupy komórek narażonych na promieniowanie jonizujące. Warszawa, Poland: Master Thesis, Faculty of Physics, Warsaw University of Technology; 2019. (supervisor: K.W. Fornalski).

24.

Krasowska

Modelowanie Monte Carlo grupy komórek poddanych działaniu promieniowania jonizującego . Warszawa, Poland: Master Thesis, Faculty of Physics, Warsaw University of Technology; 2022. (supervisor: K.W. Fornalski).

25.

Piotrowski

26.

Anandakrishnan

Varghese

Kinney

Garner

. Estimating the number of genetic mutations (hits) required for carcinogenesis based on the distribution of somatic mutations. PLoS Comput Biol. 2019;15(3):e1006881.

27.

Ponikowska

. Model Mikrodozymetryczny Komórki. Warszawa, Poland: Master Thesis, Faculty of Physics, Warsaw University of Technology; 2019. (supervisor: K.W. Fornalski).

28.

Fornalski

Dobrzyński

. Modeling of single cell cancer transformation using phase transition theory: Application of the Avrami equation. Radiat Environ Biophys. 2022;61:169-175.

29.

Hattori

Yokoya

Watanabe

. Cellular automaton-based model for radiation-induced bystander effects. BMC Syst Biol. 2015;9:1.

30.

Wysocki

Fornalski

. Computational biophysical modeling of the radiation bystander effect in irradiated cells. Radiation. 2022;2;33-51. doi:10.3390/radiation2010003.

31.

Wolff

. The adaptive response in radiobiology: evolving insights and implications. Environ Health Perspect. 1998;106(suppl 1):277-283.

32.

Feinendegen

. The role of adaptive responses following exposure to ionizing radiation. Hum Exp Toxicol. 1999;18(7).

33.

Tapio

Jacob

. Radioadaptive response revisited. Radiat Environ Biophys. 2007;46(1):1-12.

34.

Mitchel

REJ

. The dose window for radiation-induced protective adaptive responses. Dose Response. 2009;8(2): 192-208.

35.

Guéguen

Bontemps

Ebrahimian

. Adaptive responses to low doses of radiation or chemicals: their cellular and molecular mechanisms. Cell Mol Life Sci. 2019;76(7):1255-1273.

36.

Dobrzyński

Fornalski

Socol

Reszczyńska

. Modeling of irradiated cell transformation: dose- and time-dependent effects. Radiat Res. 2016;186:396-406.

37.

Dobrzyński

Fornalski

Reszczyńska

Janiak

. Modeling cell reactions to ionizing radiation: from a lesion to a cancer. Dose-Response. 2019;17(2): 1-19.

38.

Fornalski

Adamowski

Dobrzyński

Jarmakiewicz

Powojska

Reszczyńska

. The radiation adaptive response and priming dose influence: the quantification of the Raper-Yonezawa effect and its three-parameter model for postradiation DNA lesions and mutations. Radiat Environ Biophys. 2022;61: 221-239.

39.

Prise

Folkard

Michael

. A review of the bystander effect and its implications for low-dose exposure, Radiat Protect Dosim. 2003;104(4):347–355.

40.

Johansson

Carlsson

Nilsson

. Radiosensitivity of human B-lymphocytic lymphomas in vitro. International Journal of Radiation Biology and Related Studies in Physics, Chemistry and Medicine. 1982;41(4):411-420.

41.

Sanderson

BJS

Dempsey

Morley

. Mutations in human lymphocytes: Effect of X-and UV-irradiation. Mutat Res Lett. 1984;140(4):223-227.

42.

Molica

Alberti

. Prognostic value of the lymphocyte doubling time in chronic lymphocytic leukemia. Cancer. 1987;60(11):2712-2716.

43.

Seshadri

Kutlaca

Trainor

Matthews

Morley

. Mutation rate of normal and malignant human lymphocytes. Cancer Res. 1987;47(2):407-409.

44.

Prosser

Edwards

Lloyd

. The relationship between colony-forming ability and chromosomal aberrations induced in human T-lymphocytes after γ-irradiation. Int J Radiat Biol. 1990;58(2): 293-301.

45.

Macallan

Asquith

Irvine

, et al. Measurement and modeling of human T cell kinetics. Eur J Immunol. ;33(8):2316-2326.

46.

Alenzi

. The role of apoptotic proteins in patients with systemic lupus erythematosis. Egypt J Immunol. 2009;16(1):107-116.

47.

Redon

Dickey

Bonner

Sedelnikova

. γ-H2AX as a biomarker of DNA damage induced by ionizing radiation in human peripheral blood lymphocytes and artificial skin. Adv Space Res. 2009;43(8):1171-1178.

48.

Wang

Lawrence

Wan

, et al. SF3B1 and other novel cancer genes in chronic lymphocytic leukemia. N Engl J Med. 2011;365(26):2497-2506.

49.

Zheng

Cheng

, et al. c-FLIP is involved in tumor progression of peripheral T-cell lymphoma and targeted by histone deacetylase inhibitors. J Hematol Oncol. 2014 7(1):1-11.

50.

Khisroon

Gul

Khan

, et al. Comet assay based DNA evaluation of fuel filling stations and automobile workshops workers from Khyber Pakhtunkhwa province, Pakistan. J Occup Med Toxicol. 2015;10(1):1-6.

51.

Zhang

Wiencke

Koestler

Salas

Christensen

Kelsey

. Absence of an embryonic stem cell DNA methylation signature in human cancer. BMC Cancer. 2019;19(1):1-12.

52.

van Rijn

Heimans

van den Berg

van der Valk

Slotman

. Survival of human glioma cells treated with various combination of temozolomide and X-rays. Int J Radiat Oncol Biol Phys. 2000;47(3):779-784.

53.

Reszczyńska

2020. Modelowanie odpowiedzi komórkowej na małe dawki promieniowania jonizującego. PhD Thesis, National Centre for Nuclear Research. (supervisor: L. Dobrzyński; auxiliary supervisor: K.W. Fornalski).

54.

Fornalski

Reszczyńska

Dobrzyński

Wysocki

Janiak

. Possible source of the Gompertz law of proliferating cancer cells: Mechanistic modeling of tumor growth. Acta Phys Pol, A 2020;138(6):854-862.

55.

Steel

. Growth Kinetics of Tumors. Oxford: Clarendon Press; 1977. ISBN 0-19-857388-X

56.

Wheldon

. Mathematical Models in Cancer Research. Bristol: Adam Hilger; 1988. ISBN 0-85274-291-6.

57.

Dobrzyński

Fornalski

Feinendegen

. Cancer mortality among people living in areas with various levels of natural background radiation. Dose-Response. 2015a;13(3):1-10.

58.

Dobrzyński

Fornalski

Feinendegen

. The human cancer in high natural background radiation areas. Int J Low Radiat. 2015;10(2):143-154.

59.

Fornalski

KW.

Radioadaptation and radioresistance during deep space travels. The Journal of Space Safety Engineering 2022;9(3):385-389. doi:10.1016/j.jsse.2022.04.001.