Modeling Cell Reactions to Ionizing Radiation: From a Lesion to a Cancer

Abstract

This article focuses on the analytic modeling of responses of cells in the body to ionizing radiation. The related mechanisms are consecutively taken into account and discussed. A model of the dose- and time-dependent adaptive response is considered for 2 exposure categories: acute and protracted. In case of the latter exposure, we demonstrate that the response plateaus are expected under the modelling assumptions made. The expected total number of cancer cells as a function of time turns out to be perfectly described by the Gompertz function. The transition from a collection of cancer cells into a tumor is discussed at length. Special emphasis is put on the fact that characterizing the growth of a tumor (ie, the increasing mass and volume), the use of differential equations cannot properly capture the key dynamics—formation of the tumor must exhibit properties of the phase transition, including self-organization and even self-organized criticality. As an example, a manageable percolation-type phase transition approach is used to address this problem. Nevertheless, general theory of tumor emergence is difficult to work out mathematically because experimental observations are limited to the relatively large tumors. Hence, determination of the conditions around the critical point is uncertain.

Introduction

The development of cancer in the body by transition of normal cells to cancerous ones is a complicated multistep process in which many nonlinear processes play a significant role. A living cell is a very complex biophysical system. Radiation-induced adaptive response, the bystander effect, and abscopal effects at low-radiation doses and dose rates are the key processes that need to be addressed when modeling radiation carcinogenesis.

According to the classical theory, cancer is initiated by a set of mutations in certain genes in a cell. The mutations arise because of inefficient DNA lesion recognition and/or repair. Unrepaired lesions that result from the DNA replication errors and the activity of metabolic free radicals can produce spontaneous mutations occurring at the rate from about 1·× 10⁻⁷ to 5·× 10⁻⁶/gene/cell/year. ¹ This relatively broad range, which is important for further considerations, may be due to variable output and/or activity of repair proteins. Their activity, in turn, can be modified by exogenous stimuli.

Exogenic lesions, caused, for example, by exposure to intermediate and high doses of ionizing radiation, may lead to the development of cancer through radiation-induced genetic and epigenetic changes. Within this framework, biological effects of ionizing radiation are modeled by a step-by-step introduction of key processes that lead from single changes in the DNA to a full-blown cancer. In contrast to exposures at intermediate- and high-radiation doses, absorption of low-radiation doses, especially when delivered at low dose rates, is unlikely to produce multiple irreparable DNA lesions but still alerts the DNA damage surveillance system. This results in a stimulated repair of numerous DNA lesions in genes, including those associated with cell replication and metabolism. In this case, the increased repair capacity of the cells (and therefore of the whole organism in the case of whole-body exposure), which is manifested by a decreased overall rate of fixed mutations in the DNA, can translate into reduced risk of neoplastic transformation of cells and of cancer development. The degree of natural protection stimulated by low-radiation doses depends on the type of radiation, its dose, and dose rate.

Within a standard adaptive response study design, a small priming dose is used to upregulate adaptive response mechanisms (which represent a mild stress response and in vivo can involve a hierarchy of natural protective mechanisms). ^2,3 A large challenging dose is then administered (usually shortly after the priming dose), and biological end points (eg, the rates of cell deaths, mutation or neoplastic transformation) of such a combined exposure are compared to the ones when only the challenging dose is used. A reduced frequency of adverse biological effects in the presence of the priming dose indicates a rapid adaptive response (i.e., rapid adaptation) induced between the 2 exposures, which may involve epigenetic changes. ⁴ The priming dose can be brief or protracted for these effects to take place. In addition, the priming dose can lead in vivo to reduction in the rates of mutations ⁵ and neoplastic transformation ⁶ to a level below the spontaneous frequency, presumably as a result of upregulation of the body’s natural defenses. ^2,3

An important role in cancer formation may be played by close-by and distant cells in a tissue through intercellular signaling. These signals are responsible for induction of bystander (nearby cells) and abscopal (distant cells) effects. However, detailed mechanisms of these effects have not yet been fully resolved and there are inconsistencies in their understanding. For example, some authors reported that in studies of the bystander effect, signals from irradiated cells to a unirradiated cells surrounding an irradiated one exacerbate lesions in the latter. ^7,8 However, other studies demonstrate elimination of such lesions (e.g., Mothersill and Seymour ⁹ ) through the induction of apoptosis (eg, self-destruction of both hit and not hit-transformed cells). The latter is now considered a different form of adaptive response and relates to stimulated upregulation of natural protective mechanisms. ² In the case of the bystander effect, one can try to model it (eg, Hattori et al ¹⁰ and Khvostunov and Nikjoo ¹¹ ) and use Monte Carlo simulations (eg, as in the article by Fornalski et al ¹² and Hattori et al ¹⁰ ), and the result will clearly depend on the employed model. It is believed that mechanisms of bystander effects and adaptive response (commonly associated with low-dose exposures) are basic components of the cellular homeostatic response. ⁷ Adaptive response to radiation has been described theoretically and modeled in a number of ways. However, none of the currently available quantitative risk models have covered the path from the deposition of radiation energy in a cell to a developed (i.e., neoplastically transformed) cancer cell. As will be discussed in the sixth section, in order to fully describe this process, it is evident that any change from a simple collection (set) of cancerous cells to a more complex system of a malignant tumor must be associated with a basic reorganization of the initial set into a new entity, the properties of which cannot be readily and uniquely derived from the properties of the initial system. Until a more precise understanding of the underlying mechanisms is reached, we will focus only on description of hit cells. Thus, one has to accept the fact that the calculations presented in this article will have to be corrected for the 2 aforementioned effects once their mathematical form and relative strengths are worked out.

Describing (or modeling) the development of radiogenic cancer from the initial radiation-induced/exacerbated genetic and/or radiation-exacerbated genetic and/or epigenetic changes to a clinically detectable neoplasm may be regarded as an unattainable task. Moreover, since there are different ways for cancer to arise, different conceptual models have to be introduced. One of these models was introduced by Hanahan and Weinberg in their 2 seminal articles. ^13,14 Their conceptual model focused on the hallmarks of cancer which provide guidance on the key processes that should be addressed when developing a quantitative, mechanism-based model for the development (risk) of cancer.

We propose a possible biophysical interpretation of the processes of creation of radiation-induced changes in the DNA and the ensuing mutations in exposed cells. Transformation of a mutated cell into a neoplastic one is also discussed. In all these processes, the adaptive response mechanism (which has been proposed and successfully used in earlier theoretical studies) is implemented and discussed. Most of the mathematical formulae will be presented as proper probability functions that can be used in Monte Carlo simulations.

The article is organized as follows: In the second section, we address the possible outcomes of exposing a cell to ionizing radiation. In the third section, the path from lesions to mutations is described, mainly basing on the Random Coincidence Model—Radiation Adapted (RCM-RA). ¹⁵ The fourth section describes adaptive responses of cells exposed to acute or protracted (continuous) irradiation. This leads to considerations of transformation of a mutated cell into a cancer cell, including the relationships between the dose rate and the number of mutations on the number of the developed cancer cells (fifth section). A more detailed description of the inception of a tumor from precancerous cells is the subject of the sixth section. Conclusions are presented in the final section of the article. The general idea of this modeling utilized by us is presented as a flowchart in Figure 1. Every step is described in each section of the article.

Figure 1.

Flowchart of the model used. The meanings of symbols are explained in Table 1 and in the text of the article.

Table 1.

Summary of Main Probability Functions in the Presented Model.

Probability Function	Described Process or Explanation
Phit	radiation DNA hit in a cell
PL	Creation of the DNA damage from the hit
S1	initiation of single base change or a SSB per unit of time
St	any of all considered types of DNA breaks occurring in time
Pm	mutation creation per unit of time
PR	repair process, reduction of number of lesions per unit of time
PAR	radiation-activated mechanism of the adaptive response
PNR	natural repair of DNA lesions
Pc(m)	cell transformation into a cancer cell due to number of accumulated mutations in time

Creation of Lesions in a Cell After Deposition of Radiation Energy

The interaction of ionizing radiation with matter depends on the type R of radiation and on its energy, E. This interaction is typically described by the cross-section for a given process, σ (R, E). The cross-section of interaction of ionizing radiation with matter (or a cell) depends on many physical effects. As an example, an interaction of a single photon with matter is briefly described in Appendix A.

Ionization events in the cell may lead to several biophysical effects. Here, only lesions that may lead to neoplastic transformation will be considered, such as radiation-induced single-strand breaks (SSBs) and double-strand (DSBs) breaks will be considered.

As an example, for electrons of energies 5.6 eV and 9.6 eV, the values of SSB cross-sections are 2.4·× 10⁻¹⁴ cm² and 4.7·× 10⁻¹⁴ cm², and those of the DSB cross-sections are 4.1·× 10⁻¹⁵ cm² and 4.5·× 10⁻¹⁵ cm², that is, the former are by an order of magnitude smaller than the latter. ¹⁶

In addition to the cross-sections, one needs also to consider the probability of a radiation hit at the DNA, P _hit, which depends on the flux of impinging particles, and the surface density (number per cm²) of the elementary objects (e.g., the DNA or its secondary structures) to be hit. The probability of a hit can be combined with the dose D (i.e., dose rate, $\dot{D}$ , multiplied by time of exposure) absorbed by the object. Such a dependence may be linear at a low dose; however, one could postulate that it must saturate at a high dose. Thus, the total probability function describing the creation of a DNA damage (lesion) after the radiation hit, P_L , can be described as:

$P_L=A\text{σ}{P}_{\text{hit}}=A\text{σ} (1-{\text{e}}^{-cD}),$

1

where A is a normalization constant and c denotes a scaling constant. Because cross-sections are expressed in barns b (1b = 10⁻²⁴ cm²), constant A must be expressed in cm⁻². Obviously, it is the product of the thickness of the target (in cm) and the numerical density of the interacting objects, (eg, number of cells per cm³). The concept of P _hit in Equation (1) was originally used by us in the Monte Carlo chain cellular model, ¹⁷ but the validity of this formula has not been verified. Now, such a validation is presented in Appendix B.

Phenomenological Descriptions of Lesions and Mutations in Irradiated Cells

When a damage to the DNA is identified as a single lesion with probability P_L , one can consider the probability of creating a mutation resulting from an unrepaired or improperly repaired lesions. The process of mutation creation may be described by a polynomial, presumably dependent on the linear energy transfer (LET) of a given type of radiation R. ¹⁸ In particular, one can assume that the probability of a mutation caused by a mixed radiation field may have the form:

$p_1=1-\text{e}{-\sum _R^{ }\sum _{i=0}^{ }{a}_{i,R}{D}_R{ }^i}^{ },$

2

where D_R represents an absorbed dose of radiation R, and a_i,R are experimentally derived parameters for a given type of radiation R. Obviously the units of a_i , _R must ensure the dimensionless product in the exponential function. At low doses, Equation 2 can be approximated by:

$p_1=\sum _R\sum _{i=0}{a}_{i,R}{D}_R{ }^i.$

3

Equation 3 is used by International Atomic Energy Agency in biodosimetric standards, ¹⁹ where aberrations (such as dicentrics) are used to assess the dose received by the irradiated person. According to conclusions by Kellerer and Rossi ¹⁸ and Szłuińska et al, ²⁰ this probability (p₁ ) is linear for high-LET radiation, such as neutrons (i = 1), and linear-quadratic (i = 1, 2) for low-LET radiation, such as γ- or X-rays.

A different approach was proposed by Fleck et al ¹⁵ who presented a biophysical model initially considering the probabilities of SSBs and DSBs induced in the DNA either one after another (SSBs) or during a massive attack of reactive oxygen species (ROS) or ionizing radiation (DSBs). In their RCM-RA model, the dose rate serves as a crucial parameter in the role of ionizing radiation in this process. Assuming that the metabolic chemical burden production rate is C (per unit time) and the dose rate is $\dot{D}$ , the probability of a single base change or an SSB per unit of time is (original abbreviations):

$S_1= \text{α}C+\text{β}\dot{D},$

4

where α and β are weighting factors. The first term on the right side of Equation 4 describes the probability of damage (per unit time and a nucleotide) arising from the natural metabolism. The second term describes a similar effect, caused by ionizing radiation. The β constant (in Sv⁻¹ times the time unit) has also been calculated (for low-LET radiation only) by Fleck et al ¹⁵ as seen in their equation (A.10). The α coefficient is dimensionless.

Let the average time needed for error-free repair be τ. It may be expected that this must depend on the efficiency of repair enzymes. The average number of repaired lesions within τ is:

$S_2=\left(\text{α}C+\text{β}\dot{D}\right)\text{τ}.$

5

Consequently, the probability per unit time of the development of a DSB as a result of the sequential production of 2 SSBs closely related in space and time should be proportional to:

$S_t= {\left(\text{α}C+\text{β}\dot{D}\right)}^2\text{τ}.$

6

In this model, ¹⁵ the rate of the poorly repairable DSBs is linear in the time needed to repair single lesions. Even if this assumption may not necessarily hold true, it may be accepted as a first approximation. The kinetics of a DSB repair in humans was recently considered in greater detail by Jain et al ²¹ who demonstrated that irradiation at low dose rates increases the efficiency of such a repair. This finding is important in view of the level of complexity of such repair of the naturally created DSBs. This involves a synchronized action of dozens of proteins involved in the 2 repair pathways, homological and nonhomological, occurring with different accumulation speeds. In the meta-analysis carried out by Kochan et al, ²² the kinetics of the DSB proteins behaves with time t as $1-\text{exp}(-t/{\text{τ}}_1)$ , where the characteristic time ${\text{τ}}_1$ is the inverse accumulation speed, which could be modified by adaptive mechanisms.

Assuming that the production of the repair enzymes increases up to the saturation at a certain equilibrium, Fleck et al ¹⁵ argue that it is reasonable to assume that the average repair time must decrease in a manner inversely proportional to $(1+\text{δ}\dot{D})$ , where $\text{δ}$ is a coefficient related to the enzyme production rate, possibly dependent on LET. This assumption stems from the following reasoning: The higher the radiation dose rate, the greater the number of the induced repair enzymes per time unit. With the elevated efficiency of the repair enzymes, the probability of their presence in a close proximity to a damaged DNA fragment must increase. Hence, the average repair time of the damaged fragments should decrease with increasing dose rate. Equation 6 will then take the form:

$S_t= {\left(\text{α}C+\text{β}\dot{D}\right)}^2\frac{{\text{τ}}_0}{1+ \text{δ}\dot{D}},$

7

where τ₀ denotes the characteristic repair time for the nonradiation-induced lesions. At a sufficiently large αC, this equation exhibits an apparent hormetic-like dip of S_t at the low-dose rate exposure:

${\dot{D}}_{\text{min}}= \frac{\text{δα}C-2\text{β}}{\text{βδ}} .$

8

Instead of reliance on the postulated Equation 7, one could alternatively assume:

$S_t= {\left(\text{α}C+\text{β}\dot{D}\right)}^2{\text{τ}}_0{\text{e}}^{-\text{δ}\dot{D}},$

9

which, at small values of $\text{δ}\dot{D}$ , is not very different from Equation 7. Equation 9 represents a reverse situation: Instead of observing the minimum, the function in Equation 9 exhibits a maximum at the value:

${\dot{D}}_{\text{max}}= \frac{2\text{β}-\text{δα}C}{\text{βδ}} .$

10

This shows that the final description of S_t is very sensitive to assumptions, so one should be careful with postulating a definite formula for the dependence of repair time on the dose rate. In fact, the main assumption which led Fleck et al ¹⁵ to propose such a dependence (Equation 7) was that production of the repair enzymes increases linearly with dose rate, which may not necessarily be the case. To conclude, we note that Equation 8 offers a limit on the coefficient δ, namely, δ must be ≥2β/(αC) if one accepts the shortening of τ with the dose rate as in Equation 7, and <2β/(αC) if one accepts Equation 9. Apparently, since the number of DNA lesions should initially increase and decrease only after maximal accumulation of the repair enzymes, ²² Equation 9 could also be accepted based on such phenomenological considerations.

Let us note that at very low dose rates the S_t values behaves as:

$S_t\sim \text{const} (1-\text{δ}\dot{D}),$

11

where $\text{const}={(\text{α}C)}^2{\tau }_0$ . Let us also note that the positive value of ${\left(\dot{D}\right)}_{\text{min}}$ may equally well bind any of the 4 constants appearing in the Equation 8. The hormetic-like minimum observed in Figure 1 of Fleck et al ¹⁵ can be explained based on the assumption of reduction in repair time with the dose rate as in Equation 7 without considering that the repair time may vary with time after irradiation. Last but not the least, at high dose rates, S_t becomes either proportional to the dose rate, if Equation 7 is used, or tends to zero, if Equation 9 is used. In the former case, one observes a linear no-threshold (LNT) behavior, while Equation 9 demonstrates the decreasing probability of DSBs production due to the shortening of the time of the first SSB repair. In such a situation, a DSB would be expected to mainly be produced by a mechanism different than the consecutive induction of 2 juxtaposed SSBs. At very low dose rates or at small value of δ, the use of either Equation 7 or Equation 9 practically yields the same results.

Fleck et al ¹⁵ assumed that the repair of double lesions of all possible kinds cannot be successful (is error-prone) and that the rate of their appearance is proportional to the dose rate. With this assumption, Fleck et al ¹⁵ showed that the proportionality constant has the form $\left(1/f_{\text{nuc}}\right)\overline{z_F}{\text{β}}^2$ , where f _nuc denotes the average fraction of the volume of a cell nuclei (about 0.3), and $\overline{z_F}$ is the mean specific energy per event deposited in a critical volume. Thus, the modified S_t of Equation 8 is:

$S_t= {\left(\text{α}C+\text{β}\dot{D}\right)}^2{\text{τ}}_0{\text{e}}^{-\text{δ}\dot{D}}+\left(1/f_{\text{nuc}}\right)\overline{z_F}{\text{β}}^2\dot{D}.$

12

It is easy to check that the second term on right-hand side of this equation has the dimension of the inverse of time, for example, 1/s. Equation 12 can still be supplemented by a term accounting for cell killing. Again, this model leads to LNT-type relationship at high dose rates and should not be used in the considered range of dose rates. In the modification of the RCM-RA model ²³ published 2 years later after the report of Fleck et al, ¹⁵ the average deposited energy $\overline{z_F}$ was substituted by ${\overline{z_F}}^{+}= \frac{D}{1-{\text{e}}^{-D/\overline{z_F}}}$ , where D denotes the person’s lifetime dose (whether it includes natural background radiation is not clearly stated in the cited article), and ${\overline{z_F}}^{+}$ corresponds to the “mean of the specific energy deposited in the affected cell volumes, that is, volumes that have experienced at least one energy deposition event.”

Equation 12 can be modified by taking SSBs into account, using the same notion as in Equation 3 (or Equation 2 in general). However, SSBs are usually efficiently repaired, and the unrepaired DSBs dominate, so, for practical reasons, there is no need to make such a modification.

Using Equation 12 with the first term as in Equation 7, the authors of the RCM-RA model ^15,23 obtained an almost perfect fit to the Cohen data, ²⁴ corrected for smoking, which means that within the scope of the cited articles, SSBs appear to play a minor role in the creation of mutations and the ensuing neoplastic transformations of cells. The authors of the RCM-RA model got a relatively shallow minimum of the lung cancer mortality after exposures at a very low dose rate (roughly 1 mGy/y) and the apparently steady increase at higher dose rates. Such an increase has been attributed to the differential inhibition of the body’s natural anticancer defences. ⁴ Within the indicated article (which introduces a hormetic relative risk [HRR] model), stochastic thresholds for inhibition apply. In view of some other reports on the effects of radon exposures, ²⁵ this increase in the considered range of doses/dose rates is questionable.

For a more complete description of the development of a mutation, one needs to include the probability of the creation of a lesion, P_L , which is represented by Equation 1. As mentioned earlier, some additional appropriate repair mechanisms need to be accounted for in that modified formula as well to allow for a reduced probability of the creation of a mutation. Only after considering the expected reduction in the number of lesions by the repair of the DNA damage sites (such as SSBs and DSBs), one can reliably calculate the expected number of point mutations. In addition, as discussed subsequently, one mutation is not likely to produce a cancer cell. Thus, the joint probability function of mutation’s creation per unit of time can be presented as:

$P_{\text{mutation}}=P_L\cdot p_{1·}·\left(S_t-P_R\right),$

13

where p₁ represents probability of stable mutation and P_R is the general probability function (per unit of time) that describes repair mechanisms additional to the one already used in S_t (time dependent probability of creation of unrepairable DSB lesion with included repair mechanism of SSB). Let us recall that probabilities P_L and p₁ are dimensionless. Equation 13 is consistent with (but different from) the Feinendegen’s model of dual action. ^26,27 The probability of the occurrence of detrimental effects, S_t , can be described by Equation 12 or by other forms, such as the one presented by Dobrzyński et al. ²⁸

The repair probability, P_R , can be generally composed of 2 main components (P_R = P_NR + P_AR ): natural repair of the DNA lesions (cell age dependent and possibly genetically determined), P_NR , and the radiation-activated mechanism of adaptive response, P_AR . That these 2 repair mechanisms can be simply added is the assumption only, valid when the radiation-induced repair is not making use of the natural protective mechanisms. If this is not the case, P_AR should be considered as P_NR (1 + R_I ), where R_I is describing stimulation of the natural protection due to radiation. P_NR can likely be approximated by the inverted sigmoidal function (the Mehl-Avrami equation), where the repair possibilities decrease over time:

$P_{NR}=C {\text{e}}^{-a K^n},$

14

where C and a are scaling constants and K is the normalized age of the irradiated cell, that is, the real time divided by an accepted characteristic time constant. The parameter n is of crucial importance because it is determined by downregulation of the repair enzymes with the cell’s age.

A large fraction of unrepaired cells also undergo mitotic death or apoptosis (programmed cell death) and thus do not contribute to the mutation load. ²⁹ In the context of oncogenesis, both types of cell death offer a “successful” resolution of ineffective repair of the DNA damage, especially after a small number of dead cells can be tolerated. Moreover, it seems that removal of already transformed cells is also possible after irradiation at low doses through intercellular apoptotic signaling. ²⁹

Adaptive Response to Ionizing Radiation

The adaptive response phenomena include triggering of repair mechanisms, especially in the DNA, after irradiation of cells, tissues, or whole organisms. The number and efficiency of the activated repair enzymes are associated with the number of ionization events, which depend on the dose and dose rate. Adaptive response is assumed to be reasonably well accounted for using equations similar to Equations 2 or 3. The main implication of the latter is that the efficiency of repair is expected to grow continuously with the dose of radiation. However, since the effectiveness of repair saturates at high doses, this equation is expected to be reliable only for low doses.

Additionally, one can assume that the efficiency of the repair enzymes decreases exponentially with increasing dose. In what follows, we shall use an exponential decay in which:

$p_2={\text{e}}^{-bD}.$

15

Finally, the dose/dose-rate related probability of the effectiveness of the repair enzymes can be described as a product of Equations 3 and 15:

$P\left(D\right)=p_1{p}_2=\left(\sum _R\sum _{i=0}{a}_{i,R}{D}_R{ }^i\right) {\text{e}}^{-bD}.$

16

As already mentioned, all Equations 2, 3, 15, and 16 are connected with the dose/dose–rate relation.

With respect to the time dependence, one can assume, as the first approach, that the number of repair enzymes and their effectiveness increases with time (after an initiating event) with a probability:

$p_3={\text{μ}}_0+\text{μ} t,$

17

where μ describes the enzyme production rate after a pulse of radiation. This assumption should be not far from reality, especially at short times after the irradiation.

If one considers a single radiation pulse only, the effectiveness of the activated enzymes, after the initial rise in their concentration, must also decrease with time with a certain time constant (lifetime), 1/λ. ¹⁹ Were the probability of such a decrease per unit of time constant, this decrease would be described by:

$p_4={\text{e}}^{-\text{λ} t},$

18

which finally would lead to the overall time dependence:

$P\text{(}t\text{) = }{p}_3{p}_4=\left({\text{μ}}_0+\text{μ}t\right){\text{e}}^{-\text{λ}t}.$

19

The general shape of Equation 19 is similar to the shape of Equation 16 and can be generalized in an analogous manner:

$P\left(t\right)=\left(\sum _{n=0}{\text{μ}}_n{t}^n\right){\text{e}}^{-\text{λ} t},$

20

where the index n may be of noninteger type as the proportionality of p₃ with time still remaining an arbitrary assumption. For practical reasons, however, the simplified form given by Equation 19 is preferred.

Finally, the joint probability function of the adaptive response should be dependent on both the dose-rate, $\dot{D}$ , and the time, t. Obviously, the product of these 2 parameters is the absorbed dose. One should also note that at high doses the consideration of a time-dependent adaptive response makes no sense because of the smallness or non-existence of adaptedness. ³⁰ In numerical calculations, one introduces time steps, k ∈ {1,…, k _max} and the dose per unit time step (D), that is, the dose rate rather than the dose. The value of the time step has to be chosen independently. It seems convenient to use the time step equal to τ as introduced in Equation 5. As indicated earlier, both variables can be used independently in 2 different equations, depending on the context. Thus, the simplest forms of the appropriate functions are:

$P\text{(}D{\text{) = α}}_1{D}^n{\text{e}}^{-{\text{α}}_{2}D},$

21

$P\text{(}k{\text{) = α}}_4{k}^m{\text{e}}^{-{\text{α}}_{3}k}.$

22

Let us note that the normalization constants, α₁ and α₄, are dependent on the remaining parameters n and m (higher than 1 to obtain a hunchbacked shape of the curves) so that α₁ = α₁ (n, α₂), and α₄ = α₄ (m, α₃). The true dependence is determined by the assumed ranges of D and k, respectively. This approach was successfully used in the Monte Carlo modeling, where the joint probability function of the adaptive response was calculated in a discrete form ^12,28,30 :

$P_{AR}=C\sum _{k=0}^K{D}^n{\left(K-k\right)}^m{\text{e}}^{-{\text{α}}_{2}D-{\text{α}}_3\left(K-k\right)},$

23

where C represents a normalization constant and K the cell’s age given as the number of elementary time steps. This equation may be written in a continuous form ^12,28,30 :

$P_{AR}=C\underset{t=0}{\overset{T}{\int }}{\dot{D}}^n {\left(T-t\right)}^m{\text{e}}^{-{\text{α}}_2\dot{D}-{\text{α}}_3\left(T-t\right)} dt\mathrm{.}$

24

Let us note that D in Equation (23) denotes the dose per time step, whereas the dose rate in the continuous form (Equation 24) means the dose per unit of time. Obviously, such a modification requires the appropriate change in interpretation of the coefficients α₂ and α₃.

The subtle point in calculations is that one should distinguish whether the dose was delivered in a single step or continuously over a period of time. If the dose is delivered in the l’th step only, its effect at the h’th time step will be described by the simplified version of Equation 23. However, if the dose is delivered continuously from time h ₀ to time h ₁, the situation at each time step becomes more complicated and for the v’th time step:

$P\left(D,\text{v};h_0,h_1\right)=C D^n{\text{e}}^{-{\text{α}}_{2}D}\sum _{h_0}^{h_1}{\left(\text{v}-h\right)}^m{\text{e}}^{-{\text{α}}_3\left(\text{v}-h\right)},$

25

where summation runs over h and v ≥ h₁, and the dose D should be understood as a constant dose/step, that is, effectively, the dose rate. If the time step is small enough, the sum of the discrete values on the right-hand side of Equation 25 can be changed to an integral as in Equation 24:

$P\left(D,t\right)=C D^n{\text{e}}^{-{\text{α}}_{2}D} I\left(t\right),$

26

where D denotes a single-dose pulse delivered time t ago. The appropriate formulas of the function I(t) are given in Appendix C.

In experiments like those carried out by Jain et al, ²¹ the time of observation after the irradiation was close to τ or not more than a few times longer.

The abovementioned considerations are important if one wants to characterize specific situations in regions with the elevated background radiation. In the aforementioned article by Jain et al, ²¹ the level of background radiation was regarded as a priming dose relative to the additional challenge dose to the cells. The first dose was absorbed during chronic (environmental) irradiation, whereas the second dose (up to 2 Gy) was applied over 0.5 to 2 minutes, that is, in a much shorter time than the one needed for the development of any adverse reactions as well as of repair mechanisms. In a typical experiment ^31,32 demonstrating the adaptive response in cells, both priming and challenging doses were acute, that is, applied within a short period of time.

Figure 2 shows a typical priming-dose effect as a special example of the adaptive response (as modeled by us, Equations 21 and 22), when m = 1 for 2 irradiation times. For ease of the comparison, both curves were normalized to the same maximum. Figure 2 displays the case of m = 2 (Equations 21 and 22) and the irradiation applied between the 2nd and 20th time step. One can note the qualitative behavior of this response versus time which is not much different from the assumed response in each time step. During the irradiation, the response smoothly increases with every time step but does not saturate, indicating that the assumed model may not work. If it worked with chronic irradiation (eg, during environmental exposures), we should grow more resistant to it with age (the probability of adaptation saturates at older age). Obviously, our own immunological fitness deteriorates with time, so this effect must be included in such considerations. The problem is resolved when the calculated irradiation time increases. Figure 3 shows response to the irradiation time 5 times longer than in the case shown in Figure 2. The response is apparently flattening out and decreases relatively soon after discontinuation of the irradiation. Such a dependence shows that chronic exposures cause a constant adaptation of the organism to radiation which was also demonstrated in the earlier Monte Carlo studies. ³⁰ Although the strength of the maximal adaptive response is limited, it can still be substantial. Therefore, as observed by Jain et al, ²¹ inhabitants of regions with a substantially elevated background radiation can indeed present a higher radio-resistance.

Figure 2.

Normalized adaptive response (as modeled by us, Equations 21 and 22). Irradiation time 2-5 steps (red), and 2-20 steps (blue) for m = 1 (figure on the left) and m = 2 (figure on the right).

Figure 3.

Same as Figure 2 for longer irradiation time (1-100 time steps).

If the dose rate becomes too high for the enzymes to perform the error-free repair, the constant parameter related to the dose rate (e.g., b in Equation 15) should have the meaning of the inverse of characteristic dose rate which describes the effectiveness of the enzymes. In a more restrictive reasoning, one should bear in mind that the formula like Equation 2 may be different for low- and high-dose rates. We know that different groups of genes are involved in repair actions in these 2 regimes, so to stay on the safe side one has to limit our considerations to low-dose rates. Thus, the fundamental background of the adaptive response effect is described by Equations 21 and 22, with the most general form of:

$P\left(\text{ξ}\right)=a {\text{ξ}}^n{\text{e}}^{-\text{λ} \text{ξ}},$

27

where $\text{ξ}$ may denote the dose, dose rate, as well as the time. Were the validity of Equations 2 and 3 questioned, Equation 27 would still look reasonable. The hunchbacked shape of Equation 27 is commonly encountered in the literature. For example, Feinendegen found that the probability of the induction of adaptive response should be given by the probability distribution function with the maximum at low doses and the strongest effect being apparent after some period of time. ³³ The shape of this simple function is governed by 2 parameters, n and λ, only.

The maximum value of Equation 32 is attained at:

${\text{ξ}}_{\text{max}}= \frac{n}{\text{λ}},$

28

where it reads:

$P_{\text{max}}\text{ = }a{\left(\frac{n}{\text{λ} \text{e}}\right)}^n.$

29

In a special case of chronic irradiation, one can easily calculate the mean lifetime or the mean survival fraction of the repair enzymes (or their effectiveness):

$<\text{ξ}> = \frac{{\int }_0^{\infty }{\text{ξ}}^{n+1}{\text{e}}^{-\text{λξ}}d\text{ξ}}{{\int }_0^{\infty }{\text{ξ}}^n{e}^{-\text{λξ}}d\text{ξ}}.$

30

Because:

$P_{\text{total}}\text{ = }\underset{0}{\overset{\infty }{\int }}a{\text{ξ}}^n{\text{e}}^{-\text{λ ξ}}d\text{ξ}=a\frac{n!}{{\text{λ}}^{n+1}},$

31

the mean lifetime $<\text{ξ}>$ becomes equal to:

$<\text{ξ}> = \frac{n+1}{\text{λ}}$

32

One can note that the chronic low-rate irradiation can be treated as an infinite series of small radiation pulses. Indeed, integrating the sequence of Equation 27 from time zero to infinity after the irradiation time, one obtains:

$P_t=\text{α}\underset{0}{\overset{\infty }{\int }}\left(t+\text{θ}\right){\text{e}}^{-\text{λ} \left(t+\text{\hspace{0.05em}θ}\right)}d\text{θ}=\text{α}\frac{1+\text{λ}t}{{\text{λ}}^2}{\text{e}}^{-\text{λ}t}.$

33

This result shows that in spite of a linear increase in total dose at longer times, this probability exponentially decreases with time with the rate of the initial reaction to dose, that is, 1/λ. Initially, the probability increases although nonlinearly. The same effect is observed when the leading coefficient in Equation 27 is changed to t ² . Then P_t changes to:

$P_t=\text{α}\underset{0}{\overset{\infty }{\int }}{\left(t+\text{θ}\right)}^2{\text{e}}^{-\text{λ} \left(t+\text{\hspace{0.05em}θ}\right)}d\text{θ}=\text{α}\frac{2+2\text{λ}t+{\text{λ}}^2{t}^2}{{\text{λ}}^3}{\text{e}}^{-\text{λ}t}.$

34

Neoplastic Transformation of Mutated Cells

While the description presented earlier was given in terms of dose rates, the cumulative dose itself can be considered as well. In fact, so far, the only need for the time variable has been to address repair of individual SSBs and the adaptive response. This repair time is stochastic rather than deterministic and is relatively short (according to Fleck et al, ¹⁵ it takes about 40 minutes), so one speaks about very low doses when one considers low dose rates. At such doses, the epigenetic term in, for example, Equation 4 dominates, and therefore the second term on the right-hand side of Equation 11 can make a difference. This term strictly relates to a specific cell response to irradiation: production of the hard-to-repair DSBs. This response, however, even if happens in individual cells, in tissues should also strongly depend on the time elapsed since the irradiation. At a constant dose rate, the number of the repair-resistant DSBs should increase with time as should the number of the mutated cells. In Equation 11, the second term reflects the LNT approach. Thus, one must take into account that the organism counteracts a defective DSB and other lesions in tissues using repair mechanisms (natural and adaptive responses) as proposed in Equation 13.

Neglecting cooperation between cells, Fleck et al ¹⁵ suggested that the time-dependent generation of cells with the first mutation (the number of cells per person at time t which incurred 1st mutation, M₁ ) should be governed by the equation:

$\frac{d{M}_1}{dt}= \left(B_0{M}_0-B_1{M}_1\right)P_m$

35

where P_m denotes P _mutation, see Equation 13. M ₀ in this equation denotes the number of nonmutated cells, while B ₀ is interpreted as the “number of critical DNA bases in critical codons of all tumour associated genes per cell.” According to the Human Genome Project (note 1), a human genome contains about 25 000 coding genes composed of approximately 3 billion DNA base pairs. It seems that the genome includes 291 cancer-associated genes, and more than 1% of all genes are thought to be involved in carcinogenesis. ³⁴ Hence, about 1% of all the DNA bases are likely to represent such a critical value of B ₀.

For cells of the same tissue, one can assume that B ₁ (the value similar to B₀, but after the first mutation) should not be much different from B ₀. It would seem that there should be a minor error if both of these coefficients were substituted by a single one, B = B ₀ = B ₁. Let us note that in the original formulation by Fleck et al, ¹⁵ the last multiplier on the right-hand side of Equation 35 is S_t . In order to preserve our reasoning, this function was replaced by the probability of mutation, P_m , which is much closer to reality.

The solution of differential Equation 35 is:

$M_1= M_0\left(1-{\text{e}}^{-B{P}_{m}t}\right),$

36

which shows at small values of time a linear growth of M ₁ with time (as in equation B37 in Fleck et al ¹⁵ ) and at high t-values a saturation (equilibrium), hence $M_1=M_0$ .The saturation, however, most likely overestimates the number of single mutated cells, as $M_1=M_0$ means that the number of mutated cells is equal to the number of all cells.

To find the expected number of cells with 2 mutations, Fleck et al ¹⁵ consistently suggests equation similar to Equation 36:

$\frac{d{M}_2}{dt}= \left(M_1-M_2\right)P_m.$

37

The solution is:

$M_2= M_1\left(1-{\text{e}}^{-B{P}_{m}t}\right)=M_0{\left(1-{\text{e}}^{-B{P}_{m}t}\right)}^2.$

38

With the increase in time, the number of such cells must be smaller than M ₁. Equation 38 can next be easily generalized to the case of m mutations per cell (see Equation 40). It is important to note that according to this procedure, the number of mutated cells grows sigmoidally with time. This may indeed be expected as was shown in aforementioned article ²⁸ in which the sigmoidal dependence on dose resulted from overlapping number of linear dependencies.

It is not easy to calculate how the number of repair enzymes depends on time. However, one can assume that the growth should also be described by a sigmoidal function, so the postulated Equation 37 must be modified. Furthermore, since the number of mutations necessary for a neoplastic transformation of a cell is between 2 and 8, ^{35 -38} one can use a formula analogous to Equation 38 but with powers 2 to 8. It may be noted that in order to employ their modeling approach to fitting the Cohen ²⁴ data on lung cancer versus dose, Fleck et al ¹⁵ used the power m = 5, which seemed optimal.

The number of mutations in a cell, m, is critical for a possible neoplastic transformation to occur. One can assume that the probability of this transformation is 1 at, say, 10 mutations and may depend on the number of mutations in sigmoidal fashion using the Avrami-Mehl equation ²⁸ as:

$P_c\left(m\right)=1-{\text{e}}^{-0.0277m^k},$

39

which for m = 5 and k = 2 is close to 0.5 and saturates quickly to 1.0 at larger m values. Obviously, that form of the sigmoidal curve with its ad hoc assigned parameter value, as proposed here in Equation 39, does not follow from any first principles. ^{35 -37} Generalizing Equation 38 to the case of m mutations:

$M_m=M_0{\left(1-{\text{e}}^{-B P_m t}\right)}^m.$

40

Using Equation 39, one should get some estimation of the number of cells with m mutations that transform to cancer cells:

$N_{\text{canc}}\left(m,t\right)=M_0{\left(1-{\text{e}}^{-B{P}_{m}t}\right)}^m\left(1-{\text{e}}^{-0.0277m^k}\right).$

41

Equation 41 does not take into account any cooperative action within a collection of cells. It relates only to the creation of cancer cells from the mutated ones. As an example, Figure 4 shows contour plots in the coordinate system m-t for B·P_m = 0.01 and for the exponent k in the Equation 41 equal to 2 and 4.

Figure 4.

Relative number of cancer cells versus time, t, and the number of mutations per cell, m, for various critical exponents k = 2 (upper figure) and 4 (lower figure).

The development described by Equation 41 must terminate when the number of cancer cells, that is, the sum of N _canc (m, t) over m:

$N_{\text{cancer}}\left(t\right)= \sum _m{N}_{\text{canc}}\left(m,t\right),$

42

attains some critical value at which the voluminous tumor growth starts. Let us denote the time at which such a situation happens by t _cr. Figure 5 shows N _cancer (t) calculated under the assumption that the factor B is constant (independent of m) which to our understanding may be the case. Figure 5 shows that the calculated proliferation rate of cancer cells with time increases with the increasing critical index of cancer growth. This is reasonable as the increase in the critical index means that the rate of transformation to a cancer cell must rise. In all cases, the curves N _cancer (t) in Figure 5 exhibit a saturation and resemble the sigmoidal Gompertz curves. Of note, this saturated value, after summing up contributions from all values of m, can be calculated as:

Figure 5.

Number of cancer cells versus time. The curves in ascending order (from brown to blue) correspond to k = 1, 2, 3 and 4, calculated using Equation 42 and the assumption of B·P_m=0.01.

$N_{\text{cancmax}}=\underset{t\to \infty }{\text{lim}}\sum _m{N}_{\text{canc}}\left(m,t\right)=\sum _m{M}_0\left(1-{\text{e}}^{-0.0277m^k}\right)$

43

In a special case when m ≥ 4 and k ≥ 4, one can write that N _cancmax ≈ m·M ₀. The curves in Figure 5 are qualitatively similar to the ones obtained by an analytical approach of Dobrzyński et al, ²⁸ (see their Figure 3). The shape of all of the curves is virtually identical, differing only by a multiplication factor. These curves are, however, quite different from the ones obtained by Fornalski et al ¹⁷ who used Monte Carlo simulations of the cancer cells’ growth. As mentioned earlier, these curves can be perfectly described by the Gompertz curve; Figure 6 shows the fit of the Gompertz curve N _cancer (t) = 8.27844·exp[−6.51319·exp(−0.010028·t)] to the calculated points for the exemplary case of k = 4 from Figure 5. To the best of the authors knowledge, this is the first demonstration of the Gompertz curve (which traditionally describes the time of growth of cancer cells ³⁹ ) to be obtained from the combination of the probabilities and the basic biophysical properties considered in this article. (It is particularly noteworthy that the presented calculations, especially Equation 42, do not take into account the processes of cell divisions and deaths that could modify the curves N _cancer (t).)

Figure 6.

Calculated number of cancer cells versus time (open circles) from the special case of k = 4 from Figure 5, with B·P_m = 0.1. The fitted Gompertz function N_cancer(t) = 8.27844·exp[-6.51319·exp(-0.010028·t)] (solid line).

As was shown by Dobrzyński et al, ²⁸ in Equation 13, the tumor growth can be also described by the Mehl-Avrami type of equation, which is based on the nucleation and growth theory. It is important to understand that “it takes a tissue to make a tumour,” ⁴⁰ and that “the cancer induction is more a function of the tissue response and not a single cell response.” ⁴¹ Thus, consideration of what happens in specific cells rather than in a whole tissue is not sufficient. The tumor growth over time is governed by a critical index n showing the spatial type of the growth: linear, 2-dimensional (2D), or 3-dimensional (3D; n = 2, 3, and 4, respectively). Finally, the achieved tumor volume expressed in terms of the number of cells is characterized by:

$V=N_{\text{cancmax}}\left[1-{\text{e}}^{-\text{γ}{\left(t-t_{cr}\right)}^n}\right],$

44

where N _cancmax denotes the number of cancer cells in this volume, so N _cancmax corresponds to the maximum given by Equation 43 (see Figure 6). Similar reasoning was used by Laird ³⁹ who originally connected the Gompertz function with tumor growth. The coefficient γ must be proportional to the dose rate with the same power index, n, so the argument of the exponent in Equation 44 is directly connected with the dose as indicated by Dobrzyński et al. ²⁸

One can wonder when exactly begins the growth of a 3D cancer. This may be just a singular (critical) point as is common in phase transitions or in catastrophes. Transition to a self-organized state may also be considered. Whether or not this relates in any way to the self-organized criticality theory ⁴² is not yet clear.

In fact, it can be assumed that in a multistep process of carcinogenesis, ¹³ each step marks a phase transition. Consequently, the “cancer energy” landscape of a biological system can be represented by a multitude of energy valleys separated by potential barriers.

The choice of the critical index, n, is not trivial; however, it seems reasonable to limit it to 4. Because of the complex nature of a tumor, growth this index may not even be an integer. Once again, it is crucial to recognize whether one is considering the acute or chronic radiation exposure.

On Cancer Growth

Modeling carcinogenesis is an extremely difficult task because of the multitude and diversity of cancers as well as their many biological and geometrical features that have to be taken into account. As an example, one can consider a hypothetical case of a spherical tumor that obtains its nutrients from the surrounding tissues before the development of its own vascular system. The nutrients enter the tumor by diffusion, and their supply decreases with time. A solid tumor itself can contain the inner sphere of dead cells and the outer shell of live cancer cells, both quiescent and actively proliferating (eg, Aguda and Friedman ⁴³ and La Porta and Zapperi S ⁴⁴ ). In a slightly more mathematically complicated model than the ones considered heretofore by us, after neglecting the shell of the quiescent cells and assuming a constant rate of the nutrients’ consumption, one finds ⁴⁵ that the volume of the tumour changes with time according to equation:

$\dot{V}= \frac{kV}{\text{γ}}\left[1-{\left(\frac{V}{\vartheta }\right)}^{\text{γ}}\right],$

45

where γ = 2/3, k = 2a (c_z-c_a)/3, a is a scaling constant, $\vartheta =\frac{4\text{π}}{3}{\left[\frac{15}{G}\left(c_z-c_a\right)\right]}^{3/2},$ where c_z denotes the concentration of nutrients available outside of the tumor, c_a is a constant that relates to apoptotic cell death, and G denotes the rate of the consumption of nutrients. If not for the exponent γ, the Equation 45 would be identical to a logistic equation (if γ = 1) with limiting value of the volume.

A more advanced approach to spherical tumors is presented by Jiang et al. ⁴⁶ Valuable reviews of other analytical models of tumor growth can be found in the already cited book by Aguda and Friedman ⁴³ and a recent monograph by La Porta and Zapperi. ⁴⁴ Most of these models describe in mathematical terms the biology of cancer formation. In this respect, Equation 45 and the ones discussed by us in the previous sections present simplified approaches to this very complicated problem. A multiscale model of avascular tumor growth was also considered in detail in the book by Aguda and Friedman. ⁴³ By using the real data on the colorectal cancer, which has a spheroid shape, these authors showed excellent agreement of their experimental results with those of the Monte Carlo calculations shown in Figure 5 of the work presented by Jiang et al. ⁴⁶ The time dependence of the elementary volumes of clusters has been assumed to fulfill special requirements related to the capacity of cell division. As a result, the growth of the tumor volume turned out to be fairly well described, while the growth curve could be fitted with the Gompertz function. On the other hand, the solution of Equation 45 leads to time dependences with shapes similar to the logistic curve as well as to the one showed in Figure 6.

The models presented thus far may be useful in characterizing the time dependence of the carcinogenic process. They rely on the following simplistic reasoning: the DNA in cells is attacked by ionizing radiation (which is our focus) as well as by free radicals (produced during normal aerobic metabolism) that evoke lesions in the DNA structure. If unrepaired, these lesions may be passed on to the next-generation cells and give rise to mutations which, when expressed in proto-oncogenes and tumor suppressor genes, may lead to neoplastic transformations of cells.

As indicated earlier, the existence of radioadaptive responses induced by low-dose irradiations invalidates LNT model employed as a basis for radiation protection regulations. ³ This is because radiation doses used to demonstrate the adaptive response (a small dose followed by a large dose) are not additive as required by the LNT model. Moreover, according to this model, potential mutations and neoplastic transformations caused by absorption of low-radiation doses add to the number of spontaneously produced mutations and transformations. However, actual data show that exposures at small radiation doses downregulate rather than increase the amount of such spontaneous effects. ^2,3,5,47

It is now commonly accepted that cancer indeed arises from a single cell transformed through a series of genetic mutations, epigenetic events, and environmental determinants that cause and sustain ectopic expression of growth-related genes (see the reviews by Kreso and Dick ⁴⁸ and Islam et al ⁴⁹ ). The cardinal property of this single cell is its “stemness,” that is, the capacity for self-renewal and multilineage differentiation into subclones of daughter cells that, after further genetic and epigenetic changes, produce heterogenous populations of cancer cells that shape the complex ecosystem of each neoplasm (Allison and Sledge ⁵⁰ ).

Let’s look at this evolving ecosystem of cancer from another angle. It is clear that by passing from single cells to tissues and to organs, the organization of the system changes significantly. Since one deals with a complex system, ⁵¹ description of the changes should closely follow the rules of phase transitions (eg, Stanley ⁵² ) and of self-organization and complexity (eg, Heylighen ⁵³ and Kauffman ⁵⁴ ). As pointed out by Heylighen, ⁵³ “complex systems consist of many (or at least several) parts that are connected via their interactions. Their components are both distinct and connected, both autonomous and to some degree mutually dependent.” The other feature of a complex system is that its main units (in this case, cancer cells) are free in the sense that they can multiply or die. Obviously, this description fully reflects cancer development in the environment of normal cells and the extracellular matrix. Notably, however, a rather fundamental question of the nature of cancer cells has not been answered so far. According to Soto and Sonnenschein, ⁵⁵ a cancer represents a problem of tissue organization in which emergent phenomena (characteristic for complex systems) are of primary importance. Likewise, Mansury et al ⁵⁶ highlight the fact that “linear adding-up of individual cell behavior is invalid in the presence of the hypothesized nonlinear interaction among tumor cells and their environment. Moreover, nonlinearity would render it virtually impossible to predict the long-run dynamics of the system using a purely analytical approach.” Consequently, one must expect the discontinuities in a description of the transitions “from a photon/particle to a cell to a tissue and to a cancer.” In fact, even on the level of genes and cells, one observes emergent phenomena and the self-organized criticality (eg, Tsuchiya et al ⁵⁷ ). Notably, Mansury et al ⁵⁶ who claim that “malignant tumors behave as complex dynamic self-organizing and adaptive biosystems” also indicate that distinct phase-transition properties can be found in the number of cancer cell clusters and their temporal behavior versus the intrinsic capability of a single cancer cell to migrate. In this context, another important question is whether the once formed cancer cell stays as such until its death. ⁵⁸ Apparently, none of the earlier discussed analytical approaches addresses this question. One can also point out the remark by Prehn ⁵⁹ that a cancer may not be caused solely by mutations in the DNA, and a cancer cell may not stay as such forever (ie, can reverse to a normal cell state). Sotto and Sonnenschein ⁵⁵ suggest therefore that “it may be more correct to say that cancers beget mutations than it is to say that mutations beget cancers.”

In their seminal articles, Hanahan and Weinberg ^13,14 indicate that during carcinogenesis, neoplastically transformed cells acquire critical features called “the hallmarks of cancer.” These include growth factors self-sufficiency, insensitivity to antigrowth signaling, evasion of programmed cell death (apoptosis), limitless replicative potential, sustained angiogenesis, ability to invade and metastasize, genome instability and enhanced mutation rate, reprogramming of the energy metabolism, and evasion of immune destruction. Additionally, Frederica Cavallo and coworkers ⁶⁰ proposed 2 “immune hallmarks,” that is, the ability of cancer cells to thrive in a chronically inflamed environment and to suppress immune reactivity.

At the beginning of the 21st century, Schreiber and his colleagues described a process called “cancer immunoediting,” whereby the immune system, the most potent guardian against neoplasia, prevents cancer development at the early stages of carcinogenesis but also shapes (“edits”) immunogenicity of neoplastic cells and contributes to cancer development. ^{61 -63}

The cancer immunoediting process can be divided into 3 consecutive phases: (1) elimination during which incipient cancer cells are recognized by the alarmed innate immune system that triggers adaptive immune responses that specifically detect and destroy neoplastic cells; (2) equilibrium, when humoral and cellular immune mechanisms (eg, interferon-γ, interleukin-12, granulocytes, macrophages, and T and B lymphocytes) hold persisting cancer cells in check (cancer dormancy) but also shape the immune status of these cells and their environment; and (3) elimination during which the extant and “immunoedited” (ie, resistant to immune attack) cancer cells proliferate in the immunosuppressive environment facilitating cancer progression toward a full-blown, clinically detectable disease (reviewed in Janiak et al ⁶⁴ ).

From the perspective of this article, the most important is the third phase of cancer development. Apparently, during all phases of cancer immunoedition, there is a competition between stimulation and inhibition of cell proliferation and between dynamic disorder and order. Hence, the use of a deterministic approach, as we did in earlier sections and in Equation 45, cannot satisfy the needs: Our system is nonlinear and should not be described by linear equations as stated explicitly by Mansury et al. ⁵⁶ Moreover, there is little hope that a reductionists point of view will help us to comprehend how a cell is functioning until we fully understand the variety of molecular interactions in a normal cell and how cells function within normal tissue, an organ, and a cancer. As indicated by Saetzler et al, ⁶⁵ “the upward causation assumption completely neglects the contribution of the environment and of the emergent structure itself (by downward causation).” The popular somatic mutation theory (SMT) of cancer turns out to be insufficient to explain the variety of cancer behavior. A tissue organization field theory was proposed to better accommodate and explain the emerging experimental evidence related to cancer development. ⁵⁸

With regard to phase transition in a complex system, the natural question is what happens close to the point at which phase transition takes place? To what extent can such a transition be treated as an emergent phenomenon? Is the loss of control over tumor growth a sort of catastrophic event, such as an avalanche, ⁴² treated as an indicator of self-organization at the phase transition, or is it just a phase transition of the first or second order similar to the reentrant and other transitions of frustrated systems from disordered to ordered states. ⁶⁶ In the case of self-organization, an essential difference between what is going on before and after the transition relates to the fact that before the transition every cell more or less individually interacts with its closest tissue constituents. In contrast, after the transition, all cells work together: What happens in one place of the tumor has a direct influence on what happens at any other place. How the situation in one place will change the situation in another place of the tumor is hard to predict (possibly abscopal effects may be involved). At the beginning, all incipient tumor cells are fed by diffusion from the surrounding tissues and can proliferate. With time, however, the inner core of a tumor is formed. Probably, Equation 45 may roughly describe the evolution of a tumor volume, mainly of its outer shell. From the organization point of view, the larger the tumor, the more external, unbounded cancer cells can be accommodated on its surface, and the tumor would exponentially grow up to infinity. In any case, within the scope of the theory of complexity, one has to admit that the process of tumor growth cannot be reduced to individual interactions between cells as was possible during the equilibrium phase of the immunoediting process. Once the tumor is formed after passing the critical point, the tumor cells lose their individuality and become totally subordinated to the properties of this new entity. Of course, any stress, such as exposure to ionizing radiation, may change this self-organized behavior. Because of such complexities, one has to accept that the description of an organism cannot be reduced to interactions between its principal entities and that the organisms are subject to rules of organization and its variations as discussed by Mossio et al. ⁶⁷ From a purely physical point of view, an organism is an open system capable of exchanging energy and matter with its environment. These important problems are, however, beyond the scope of the present article.

Phase transition as discussed here may also be well understood based on the so-called percolation type of phase transition. In this case, the main assumption is that a single cancer cell is not a cancer itself, and only a group of these cells may constitute a tumor. So, at first, individual cancer cells may occasionally form contacts (or links) between each other, and the functional links may lead to the formation of clusters. Within the scope of the continuous time branching processes theory, one can calculate a cumulative distribution of the cancer colony sizes versus the number of cells in these colonies (see Figure 3.3 in La Porta and Zappero ⁴⁴ ). This distribution depends on the time of observation, qualitatively is not very different from the logistic curve, and pretty well describes the observations.

In a further development of such an intertwined network, all the clusters may fill up the space in such a way that nutrients provided to one of them can be transferred to any other—the percolation transition is achieved. Likewise, a disturbance occurring in one place can be propagated to any other place. In terms of the second-order phase transitions, one should talk about spatial and temporal fluctuations that grow, on average, below the transition point and, upon passing this point, freeze to a single ordered phase.

As mentioned in the previous section, a single cancer cell created during the neoplastic transformation of a mutated cell is not yet a cancer. Equations 44 and 45 described the tumor growth with time. To make a group of cancer cells a tumor, a certain number of them, h:

$h=\frac{N_{\text{canc}}}{N_{\text{canc}}+N_{\text{non\_canc}}},$

46

where N _{non_canc} is a sum of damaged, and mutated cells (N _{non_canc} = N _lesion + N _mut) must interact with each other and start to proliferate in a coherent way after passing a certain critical value, h_c .

With time, the value of h changes, but as long as cancer cells or their clusters are disconnected, the tumor has not yet emerged as a separate entity. As already indicated, such an emergence can be identified with a phase transition similar to the re-entrant phase transitions known from the physics of magnetics or the percolation theory. In both cases, the final formation of a given object (eg, a tumor) appears when h exceeds some critical value (eg, the number of cancer cells, h_c ) as mentioned earlier. Within the scope of the percolation theory, which refrains from purely physical or biological parameters, one needs to define the parameter that should control the occupancy of pixels (voxels) into which a given space is subdivided. This parameter should reflect not the relative number of cancer cells as in Equation 46 but rather the probability of the creation of a cancer cell at a given pixel (voxel). Let this probability be denoted by p. With time, cancer cells aggregate and form clusters that combine with each other and finally a critical state is attained: The “infinite cluster” (a tumor) is formed in which any information sent from one location in the object can reach all other locations with a consequence to the whole object—the tumor starts to behave as an entity whose behavior cannot be derived from individual properties of cells and their interactions. This is similar to the sand-pile experiment ⁴² : Although the grains of sand drop to the sand’s cone from the top onto a single point, at a certain height and radius of the cone avalanches appear in an unpredictable manner (the so-called self-organized criticality). In the case of percolation, the situation may be visualized by imagining a number of pixels or voxels into which one drops small grains or small spheres. Next, let us connect randomly any 2 such small spheres. If we repeat this procedure, the number of the connected spheres will increase, and the number of locally connected spheres (clusters) will increase. At a certain moment, the size of those clusters will start to suddenly rise. This illustrates a case of self-organization.

In a typical simulation of a percolation phenomenon, ⁶⁶ one has to choose the size of the object divided into pixels and by the “infinite cluster” one understands the cluster extending from one edge of the object to another one. Ideally, the object should have infinite dimensions so that the meaning of “infinite cluster” is literal. In our case, the situation is somewhat different. An organ is close to be infinite with respect to the size of the cells. The size of tumors in it may not be as large as the organ, for example, the lack of nutrients needed the cancer to grow. ³⁹ Nevertheless, one can still treat the maximum size of the tumor as roughly equivalent to an infinite lattice of cells, and the percolation theory ⁶⁸ with its purely geometrical statistical ingredients can be useful for the description of at least a region close to the phase transition point p_c . In the case of a site percolation on the square lattice p_c ≈ 0.593, the “infinite cluster” may exist only above p_c . Inside the organ, which represents a truly infinite lattice, one can imagine formation of more than 1 “infinite cluster,” that is, of more than 1 tumor.

According to this theory, the probability that a given cell belongs to the collection of tumor cells grows at p > p_c in a critical way, and the percolation probability, P _max, that is, the fraction of the occupied sites belonging to the “infinite cluster,” is ruled by the critical index β which may in general be of the fractal type:

$P_{\text{max}}\sim {\left(p-p_c\right)}^{\beta }.$

47

This behavior is displayed in Figure 6. It is neither sigmoidal type nor logistic type; however, one should remember that the relation (Equation 47) must describe the behavior mainly in the critical region, that is, relatively close to p_c . Alas, the width of the critical region is difficult to predict.

A very illustrative example of this kind of behavior can be found in many disordered magnetic systems, such as a diluted ferromagnet which below p_c becomes a paramagnet and becomes ferromagnetic above p_c . In addition, over a certain concentration range of magnetic species (Co, Fe, and so on), the spin-glass phase can be formed. In such a system, the control parameter is temperature in Figure 7 ^69,70 —with increasing temperature, the value of the percolation probability p_c increases as well. As an example, the onset of ferromagnetism in a percolating 3D network of the random face-centered cubic alloys can be satisfactorily described within the framework of the percolation theory, where p_c ranges from 0.16 to 0.20. ⁷¹ Often, the phase diagram of such diluted ferromagnets can be described within the framework of the so-called Ising model. In the case of cancer, one encounters a more complex situation because the number of states representing cells with different degrees of lesions and mutations is much higher than the number of up- and down-spins in Ising model. This makes statistical description of the state below p_c more difficult.

Figure 7.

Schematic phase diagram of a diluted ferromagnet (eg, Au-Fe, Co-Cu, and so on) showing phase transitions between paramagnetic (PM), ferromagnetic (FM), and spin-glass (SG) phases depending on temperature and concentration X of the magnetic sample.

It should be mentioned that the sharp transition shown in Figure 8 may be smeared if one finds correlations between the cells below p_c . If such correlation exists between the states of 2 cells separated from each other by a distance R, one can find a function describing such a correlation. In a typical magnetic system, it would be described by a function decreasing exponentially with the distance. This would lead to a substantial change in Figure 8: the whole curve becoming sigmoidal-like, such as the one in figure 3.9 of Binder and Kob ⁶⁶ .

Figure 8.

Transition to percolation after passing the p _c point assumed to be p _c = 0.5 at different critical indices β.

Just below the transition point, the average distance between the cells within a cluster (correlation length), ξ, behaves as:

$\text{ξ}\sim {\left(p_c-p\right)}^{-\nu },$

48

where ν denotes another critical exponent. The values of both exponents depend on grid and slightly differ depending on case (note 2). For example, for a 3D net β = 0.418 ± 0.001 and ν = 0.875 ± 0.008. The mean cluster size also exhibits nonanalytical behavior:

$S\sim {\left(p_c-p\right)}^{-\text{γ}},$

49

with γ = 1.793 ± 0.003. In the case of 2D growth, the critical exponents change to 5/36, 4/3, and 43/18 for β, ν, and γ, respectively. Let us note that the critical exponents depend on the dimensionality of the problem and not on the microscopic properties of the system. The whole situation resembles the behavior of magnetic systems: Our p plays a role similar to temperature, percolation probability P _max—magnetization, ξ —correlation length, and S—magnetic susceptibility. Below the transition point, one can also make intuitive use of the so-called mean field approach. Namely, let the number of cancer cells be N _canc, whereas N _canc,₀ is the number of cancer cells whose growth seem to be inhibited by immunological forces and represent dormant cells. The development of cancer cells may depend on their number that form a special field due to interaction of these cells with their environment (in fact, it is the question of competing forces between natural expansion of cancer cells and the immune restrain of such a development, see Janiak et al ⁶⁴ ). This field describes something analogous to a promotion of the development of cancer with a multiplication factor, say, χ. Then one can write a simple equation:

$N_{\text{canc}}=N_{\text{canc},0}+ \text{χ}{N}_{\text{canc}}.$

50

Thus:

$N_{\text{canc}}= \frac{N_{\text{canc},0}}{1-\text{χ}},$

51

which shows how the number of cancer cells can be strongly enhanced by the interaction of these cancer cells with their environment. In an extreme case, the number of cancer cells can reach the total number of cells if χ = 1, that is, when immunological protection against cancer development breaks down. Equation 51 is a typical result obtained within the framework of the mean field theory of phase transitions. The problem to be solved is the description of the form of χ. It is not a single number but a function depending on the processes of immunological protection against cancer. In fact, all the aforementioned steps in cell behavior represent some phase transitions, although an appropriate mathematical description of these transitions is rather difficult. Consequently, one needs to define the so-called order parameter that changes upon transition. This parameter must have defined dimensionality, while one should also define whether the main interaction mechanism with environment is short- or long ranged. Besides, it seems that one deals with 3 space–time scales: the one of signal transduction, the short-range adaptive response and bystander effect, and the long-ranged abscopal effect.

Presently, the above-described critical phenomena are difficult to observe at present, as the minimum size of the detectable tumor is of the order of a few millimeters when the tumor is already formed. In spite of this obstacle, one should understand that the general properties of phase transitions, including self-organization and/or self-organized criticality, have to be included in a rigorous description of the tumor development.

Conclusions

There are several problems covered in this article. Its intention was to separate the evolution of individual cells from that taking place in a tissue. With this aim in mind, the approach of Fleck et al ¹⁵ was initially examined. Although the ideas behind their model are quite different from the ones presented previously by us, ²⁸ the final result for the number of cancer cells versus time has turned out to be qualitatively similar. Notably, in contrast to Fleck et al ¹⁵ who used dose rate, Dobrzyński et al ²⁸ employed dose in their calculations. Moreover, a more fundamental difference in the way a cell repairs its DNA lesions is presented in this article. In the model of Fleck et al, ¹⁵ the dependence of the time of repair of a single DNA lesion on the dose rate plays an essential role. Apart from the fact that a postulated modification of the repair time with increasing dose rate may be described by a different function, the model of Fleck et al ¹⁵ explains Cohen ²⁴ data on lung cancer mortality versus radon-specific activity. However, it seems that one cannot go too far with a function fitted to data, as the mortality calculated in this way soon exceeds the hormetic minimum. Such an increase is not confirmed by the numerous other data collected for much higher radon concentrations. ^25,72

Our attention has also been paid to the dependence of adaptive response on the time and the dose associated with acute or protracted radiation exposures. We demonstrate that in the case of a protracted exposure, the organism attains a certain saturation in its ability to repair lesions (Figure 5). This observation permits us to treat absorbed dose as a priming dose that allows to better toleration of higher and more challenging doses.

Increase in the number of cancer cells in an organism depends on the timing and the dynamics of the critical number of mutations in a cell needed for its transformation into a cancer cell. As demonstrated by us, this may be described by a pretty complicated function (Equation 41). However, the increase in the total number of tumor cells (Figure 5) when graphed resembles a sigmoidal shape of the Gompertz function. As previously indicated by Dobrzyński et al, ²⁸ the 3D growth of the tumor volume may generally be described by a sigmoidal curve given by the nucleation and growth according to the Avrami-Mehl theory as discussed in Dobrzyński et al. ²⁸

Finally, we have discussed at length the problem of tumor growth. Description of this process is inherently difficult. First, it essentially deals with the phase-transition phenomenon, and the type of this transition is not clear. Essentially, there is a transition from a disordered phase (represented by individual cancer cells and their clusters) to an ordered one (the developed tumor) which exhibits properties not directly connected to individual processes occurring in its basic units (cells). If this is the case, one is dealing with an emergent phenomenon, the self-organization, and the possible self-organized criticality. However, even if the phase transition can be viewed as a continuous one (of second order), it is difficult to specify its most important parameters: the order parameter and its dimension, as well as the interaction range between cells and the control parameter (usually temperature in a magnetic system or in, e.g., liquid–gas transition). Depending on these parameters, one could specify the values of critical exponents characterizing the transition. An important property of such critical indices is their universality—they are not directly related to microscopic interactions between the basic units of the system. However, when a mathematical description is given, its experimental verification can be difficult. Today, a tumor cannot be detected if its diameter is around 2 to 3 mm. This means that it is already far beyond the phase-transition region where critical properties can be observed. To illustrate the problem, phase transitions of the percolation type were considered. Additionally, a simple application of the mean field theory shows that nonanalytical behavior is to be expected if a phase transition takes place. We plan to address these problems in a future publication.