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

The priming dose effect, called also the Raper–Yonezawa effect or simply the Yonezawa effect, is a special case of the radiation adaptive response phenomenon (radioadaptation), which refers to: (a) faster repair of direct DNA lesions (damage), and (b) DNA mutation frequency reduction after irradiation, by applying a small priming (conditioning) dose prior to the high detrimental (challenging) one. This effect is observed in many (but not all) radiobiological experiments which present the reduction of lesion, mutation or even mortality frequency of the irradiated cells or species. Additionally, the multi-parameter model created by Dr. Yonezawa and collaborators tried to explain it theoretically based on experimental data on the mortality of mice with chronic internal irradiation. The presented paper proposes a new theoretical approach to understanding and explaining the priming dose effect: it starts from the radiation adaptive response theory and moves to the three-parameter model, separately for two previously mentioned situations: creation of fast (lesions) and delayed damage (mutations). The proposed biophysical model was applied to experimental data—lesions in human lymphocytes and chromosomal inversions in mice—and was shown to be able to predict the Yonezawa effect for future investigations. It was also found that the strongest radioadaptation is correlated with the weakest cellular radiosensitivity. Additional discussions were focussed on more general situations where many small priming doses are used.


Introduction
The radiation adaptive response effect (called also a radioadaptation) is a biophysical phenomenon, which may occur in the organism irradiated to low doses of ionizing radiation. This effect connects the irradiation process with the adaptation to radiation, which causes faster DNA lesion repair and reduction of DNA mutation creation (Olivieri et al. 1984;Azzam et al. 1994;Wolff 1998;Feinendegen 1999;Tapio and Jacob 2007;Dimova et al. 2008;Mitchel 2009;Guéguen et al. 2019). However, because not all irradiation conditions induce radioadaptation (Mortazavi et al. 2003;Wójcik et al. 1996), and the effect shows strong individual variability (Bosi and Olivieri 1989), it can be very difficult to predict. Therefore, many scientific investigations in radiobiology and radiation biophysics aimed at understanding this phenomenon are still going on.
Till recently, it has been assumed that the radiation adaptive response effect occurred when enhanced repair mechanisms started after the creation of damage within the DNA chain (Mitchel 2009). Nowadays, new findings show that small radiation dose(s) can affect the expression of gene transcription (Sokolov and Neumann 2015). Exposure in proper conditions can produce "an alert, triggering and altering cellular responses to defend against subsequent high dose-induced damages, and accelerating the cell repair process. Moreover, the p53 signalling pathway was found to play critical roles in regulating DNA damage responses" (Hou et al. 2015) because "the p53 pathway is composed of a network of genes and their products that are targeted to respond to a variety of intrinsic and extrinsic stress signals that impact upon cellular homeostatic mechanisms that monitor DNA replication, chromosome segregation and cell division" (Harris and Levine 2005;Vogelstein et al. 2000). Therefore, the adaptive response mechanism may increase DNA repair up-regulation, leading to the reduction of postradiation mutation frequency (Dimova et al. 2008).
There are two main processes of triggering the radiation adaptive response: via chronic irradiation or pulse single doses. In the first case, the effect can be observed e.g. within some people living in areas with a high background radiation (Scott et al. 2009;Dobrzyński et al. 2015aDobrzyński et al. , 2015b. The second case is equivalent to the priming dose effect, also called the Yonezawa phenomenon (originally for mice) (Yonezawa et al. 1990(Yonezawa et al. , 1996Matsumoto et al. 2004;Tapio and Jacob 2007), where the small priming dose (D 1 ) received prior to the high challenging dose (D 2 ) can reduce detrimental effects of the latter Wang et al. 2013;Toossi et al. 2016;Hauptmann et al. 2016). The exemplary scheme of this effect is presented in Fig. 1, where the difference between the D 1 + D 2 scheme (with time Δt between doses) and the single D 2 one is denoted as delta (δ) parameter, which is the main quantification of the generalized Yonezawa effect. One has to note, however, that the first observation of this effect (without deeper explanation) was conducted by Prof. Raper in 1940s during the Manhattan Project (Raper 1947;Cronkite et al. 1950).
The first theoretical explanation of this phenomenon was given by Dr Yonezawa and collaborators, who created the phenomenological multi-parameter model describing the effect Yonezawa 2003, 2004). The approach proposed by Smirnova and Yonezawa is a deterministic model, describing the effect of ionizing radiation on the cells of four major hematopoietic lines: thrombocytes, lymphocytes, erythrocytes, and granulocythes. It divides cells into categories depending on the stages of development (bone marrow precursor cells, nondividing maturing cells, and mature blood cells) as well as on the level of damage (undamaged cells, damaged cells, and heavily damaged cells), and describes the relationships between those categories. To describe the damaging effect of ionizing radiation on the cells, the one-target-one-hit theory is used. The model consists of a system of differential equations with many unknown (free) parameters that must be determined based on experimental data. It allows the simulation of the effect of either chronic irradiation or pulse single doses. It should be noted that this model has limited application to certain types of cells (namely thrombocytes, lymphocytes, erythrocytes, and granulocythes only) and that the large number of equations and unknown (free) parameters make it hard to apply in practice.
The presented paper proposes a novel theoretical approach to the radiation adaptive response phenomenon, leading to the Yonezawa effect. This results in a threeparameter model, which works well after parametric quantification based on experimental radiobiological data.
The model is divided into two parts: in the first part the calculations focus on the reduction of DNA damage (lesions) in quite a short time after the irradiation; in the second part late effects, namely stable mutations, are considered. The latter case is much easier to consider in experimental radiobiological research.

Adaptive response phenomenon
The proposed theoretical approach to the radiation adaptive response effect was originally presented a few years ago (Fornalski 2014) but its detailed biophysical explanation has been published very recently (Dobrzyński et al. 2019). Therefore, only a short and general mechanistic description will be presented here.
Generally, linear and non-linear effects can cause a response described by the specific hunchbacked shape (Feinendegen 2005(Feinendegen , 2016 of time-and dose-related Fig. 1 The scheme of the Yonezawa effect (also called the priming dose effect or the Raper-Yonezawa effect): the single small dose (D 1 ) generates much less mutations than the single high dose (D 2 ). However, when D 2 follows D 1 with some time distance between them (Δt), the mutation (or lesion) frequency for that D 1 + D 2 total dose is lower than for single D 2 by the exemplary value of δ = 0.73. This exemplary result was obtained using the input data: D 1 = 1 UAD (Unit of Absorbed Dose), D 2 = 5 UAD, Δt = 3 UT (Unit of Time), The parameter δ is therefore showing the percentage difference between the number of mutations (or lesions) generated by the single dose D 2 (without the priming dose) and the combination of D 1 + D 2 probability functions of the radioadaptation appearance, respectively: where D and t denotes radiation dose and time after irradiation, respectively, and {α}, {β}, ν, η are free parameters. Both functions can be merged into a single, time-and dosedependent, probability distribution of the adaptive response appearance after the single dose pulse D received t time ago (Fornalski 2014;Dobrzyński et al. 2019): where all parameters are positive ones while ν and η parameters were assumed to be equal to 2, which follows the result of pooled simplified quantification by Feinendegen (2016). Equation (2) represents the simplest version of the adaptive response per unit of time. Please note that the function p AR reaches its maximum value for the dose D max = 2/α 1 and for the time t max = 2/α 2 after irradiation. In general, Eq. (2) can be rewritten in the continuous form as (Fornalski 2014): where Ḋ corresponds to the time-dependent dose-rate and τ is the cell's (or the organism's) age. In practice, however, Eq. (3) is usually used for chronic irradiation, where Ḋ is given by some dose distribution, Ḋ(t) (Dobrzyński et al. 2019). But in the biophysical calculations it is easier to apply Eq. (2) (or its combinations) (Dobrzyński et al. 2016;Fornalski et al. 2017;Fornalski 2019). When the dose distribution is discrete, namely single dose pulses D k in time steps k, the dedicated form of Eq. (3) can be used e.g. in Monte Carlo simulations (Fornalski 2014;Loan et al. 2019). Equation (4) shall be understood as: each dose D k received K-k steps ago (K is the age) generates a single signal given by Eq. (2) extended over time, which is additive to the rest K − 1 signals described by an identical mathematical form (Fornalski 2014). The summation presented in Eq. (4) allows one to freely add every single adaptive response signal generated by each dose D k under two conditions: • the dose D k is a short pulse with the duration of t Dk ⟶ 0; • the time interval between two consecutive doses is large enough, Δt >> t Dk . (1) If some of the conditions presented above are not fulfilled, Eq. (3) shall be used instead. In particular, for two separate dose pulses, namely D 1 and D 2 , which are received with the time interval Δt > 0, one can simply denote the summarized probability distribution of the adaptive response as p AR (D 1 ,D 2 ,t) = p AR (D 1 ,t 0 ) + p AR (D 2 , t 0 + Δt). When both doses are received in the same time (Δt = 0), one shall write p AR (D 1 ,D 2 ,t) = p AR (D 1 + D 2 ,t) because they can be simply treated as the single dose (D 1 + D 2 ).

Kinetics of DNA lesions repair
Let us assume that the dose pulse, D, generates some number of fast and direct damage, N, called DNA lesions. Regarding the actual experimental findings (Rothkamm et al. 2007;Manning et al. 2013), this dependence is assumed to be linear with additional background metabolic lesions (for a zerodose situation), therefore In other words, N from Eq. (5) represents the immediate number of physical damage (lesions) generated linearly by the dose D. All parameters {μ} are calibration parameters of Eq. (5), of which values are given in the literature e.g. in direct initial lesions testing (Ward 1995) and their background level. Especially, the parameter μ 0 corresponds to the frequency of spontaneous lesions (without radiation), and µ 1 is a linear slope. One has to note, however, that the free parameter µ 0 is rather small (µ 0 << µ 1 ) and can be simply neglected. Exemplary values of {μ} parameters for damage in human lymphocytes after neutron irradiation equal µ 0 = 0.0005 and µ 1 = 0.832 (IAEA 2011;Słonecka et al. 2018).
The repair mechanisms start just after lesion (N) appearance, however, the probability of the adaptive response is strictly correlated with the number of repaired damage after a period of time, which can be presented as (Foray et al. 2005) where p AR is given by Eq. (2) for single dose pulse irradiation. It is assumed that this mechanism is the only method of repair after postradiation lesion appearance. Please note that Eq. (6) is the basis for the general function of the remaining number of lesions in an actual moment in time, N(T). The detailed calculations are presented in the Appendix 1.
Let us assume that N 0,1 denotes the initial number of damage (lesions) induced by dose D 1 in moment zero (assumed that N 0,1 = μ 0 + μ 1 D 1 , see Eq. (5)). Moreover, dose D 2 generates an additional number of damage, N 0,2 = μ 1 D 2 , in 1 3 moment Δt (where N 0,2 > N 0,1 ). 1 Figure 2 presents the time related probability of the adaptive response (p AR ) and the number of unrepaired damage in time N(T) for a single dose D 2 and the combination of priming and challenging doses, D 1 + D 2 . Please note, that the N(T) relationship is a strongly decreasing function, however, it never goes to zero due to the hunchbacked shape of the adaptive response probability function: within our model, the stronger the repair processes are, the closer to zero the N(T) function is at a large T. This is consistent with experimental findings, see Fig. 3 based on the paper by Müller et al. (2001).
The Yonezawa (priming dose) experimental scheme, which is presented in Fig. 1, can be quantified as where N 2 corresponds to the general number of lesions (or mutations) in a single D 2 scenario, while N 1+2 is the general (2) and Eqs. (21)-(22), respectively. Plot a presents a scenario where a single (reference) dose D 2 is applied; one can observe the hunchbacked shape of the p AR function and the decrease in the number of lesions, N(t). Plot b presents a scenario with a combination of doses D 1 + D 2 , where the priming dose D 1 is given prior to the challenging dose D 2 (the latter received in the same moment as in the plot a)); one can observe that the priming dose increased the probability of successful repair (enhanced value of p AR ) which causes a stronger decrease in the number of lesions, N(t). All input data were the same as in Fig. 1. The variable t is the global time where t = 0 corresponds to the moment of D 2 Fig. 3 The residual DNA damage (lesions) represented by doublestrand breaks (DSB) in human lymphocytes related to time, after X-ray irradiation of 2 Gy. Four different relationships represents four different sensitivities to radiation, from hyper-radiosensitivity (upper grey curve) to radioresistance (lower black curve)-this last case represents the strongest adaptive response effect (Fornalski 2019). The figure was created based on the paper by Müller et al. (2001) and presentation of Feinendegen (2012) 1 Please note, that N 0,2 relationship intentionally omits μ 0 parameter to avoid double counting of background lesions/mutations in Eq. (7). number of lesions (or mutations) in a D 1 + D 2 scenario. The delta parameter, δ, is the percentage difference between N 2 and N 1+2 (Fig. 1).
After some calculations, which are expressed in detail in the Appendix 1, Eq. (7) can be written as where f are functions related to the adaptive response, see Appendix 1. Detailed forms of f functions are presented in Table 1.
The approach presented in Eqs. (8-9) describes the number of actual damage (lesions), still unrepaired, in the moment of T. This gives no information about mutations, which are created much later (>> T). One has to note, that Eq. (9) actually contains three free parameters (α 0 , α 1 , and α 2 ) plus two time variables (Δt and T, where 0 < Δt < T), which are conditions of the dedicated experiment. In practice, however, the presence of the Yonezawa effect related to fast damage (lesions) is more difficult from a radiobiological point of view, because Eqs. (8-9) vary with time. Therefore, a more popular and representative situation is described by the second part (Chapter 4) of the proposed model, namely the Yonezawa effect on mutation frequency (late effects).
The parameter δ can vary from 1 (ideal protection and full reduction of all mutations) to some minimal value δ min which can even be lower than zero (when the number of mutations after the D 1 + D 2 scenario is a little bit higher than the number of mutations after the single D 2 scenario, but still lower than the sum of independent D 1 and D 2 ). The latter case can be therefore treated as a special case of the Yonezawa effect, however, in most studies it is assumed that δ min = 0 (the case when the number of mutations after D 1 + D 2 scenario is the same as the number of mutations after single dose D 2 ), see Appendix 1 for more information.

Mutations in DNA
Here one considers the phenomena of mutations, which are permanent changes in DNA caused by mis-repaired lesions: all repair processes are finished now and all N(T) functions stabilized. That way one can consider Eq. (9), but for infinite time, T ⟶ ∞. This assumption means that the repair time is long enough to repair all possible lesions and only mutations are left. Thus, the analogous form of Eq. (9) becomes which describes the Yonezawa effect for mutations, namely late effects. Appendix 1 contains details of the above calculations. Both f functions are presented in Table 1. Exemplary results of the parameter δ for the conditions from Fig. 2 are: δ = 0.725 for lesions (Eq. (9)) and δ = 0.728 for mutations (Eq. (10)), respectively. This latter case is presented in Fig. 1. One can generally see that the difference between lesions and mutations analyses is, within the scope of the model, generally marginal for a large T (> Δt).

Two priming doses
The single priming dose, D 1 , is just a special case. Another possible case can involve two similar priming doses: the first one D 1 , the second one (D 2 ) received Δt 1 time after the first one, and the challenging dose (D 3 ) received Δt 2 time after D 2 . One can note that D 1 ≈ D 2 < D 3 . It follows from Eq. (9) and Fig. 2b that the adaptive response of the second dose (similar to the first one) results in an apparent decrease of the N(T) function as compared to the one resulting from the D 2 dose alone. Thus, one should expect that the challenging dose D 3 given at a later time starts its repair functions from smaller number of lesions than when applied at Δt after the first priming dose. In another words: two consecutive priming doses, given in proper time and value, can boost the repair process triggered by D 3 . This approach needs more laborious calculations, which are presented in Appendix 1 as well. Finally, the parameter δ can be now defined as 1 − N 1+2+3 /N 3 , so which describes the Yonezawa effect when the challenging dose D 3 is applied after two priming doses D 1 and D 2 . Additionally, f functions are presented in Table 1.
Exemplary results of calculations, based on conditions from Fig. 2 with second priming dose, D 2 = D 1 , are presented in Fig. 4.
Please also note, that Appendix 1 contains analogical calculations for multiple priming dose scenarios.

Numerical methods
Both approaches, namely Eq. (9) for lesions and Eq. (10) for mutations, can be applied to experimental data to estimate the values of α 0 , α 1 and α 2 parameters. However, the differences in results obtained by those two equations are usually marginal for a large T, as discussed previously.
The application of the model, irrespective of the equation used, needs advanced numerical methods of non-linear optimization to estimate all input parameters based on experimental output. Here two of them were used: Simplified Genetic Algorithm (SGA) as well as Bounded Limitedmemory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS-B) one. Both methods and their applications are described in Appendix 2 in details.

Lesions
Equation (9) was applied to real experimental data where human lymphocytes were irradiated in vitro (X-ray) in three consecutive radiobiological researches conducted by Fig. 4 The case from the Fig. 2 with additional priming dose D 2 = D 1 , received Δt 1 = 3 UT after D 1 . The challenging dose D 3 was received Δt 2 = 2 UT after D 2 Shadley and Collaborators Shadley and Dai 1992). These studies were selected for analysis because of a large amount of data representing the Yonezawa effect. The first experiment ) studied the effect of 3-aminobenzamide (3AB) on the repair functions present during the adaptive response; additionally, the experiment studied the effect of the conditioning dose level on the priming dose effect. The lymphocytes for testing were taken from the peripheral blood of healthy 30-40 y.o. females. The blood was exposed to low doses at 32-34 h after stimulation, the challenging dose was given at 48 h. For one group, 3AB was applied immediately after the high dose. The results showed that application of 3AB immediately after the higher dose did reverse/eliminate the repairs of the adaptive response. Based on the results, low doses of 5, 10, 50, 100 mGy are capable of inducing an adaptive response when a higher dose is applied, while doses above 200 mGy do not have the ability to adapt human lymphocytes to ionizing radiation. This is perfectly illustrated in Fig. 5, where the priming dose D 1 from the approximated range 10-100 mGy gives the highest value of δ.
The second analysed study  considered the effect of the cell cycle phase during which certain doses are applied as well as how long the effect of a priming dose can last. Low doses were applied either before 2% phytohemagglutinin (PHA) addition (G 0 ) or at times corresponding to G 1 , S or G 2 phase of the first cell cycle. Higher doses are applied either in the same or the next cell cycle (40, 48, 66, 90, or 114 h after PHA addition).
The collected data demonstrated how the response to lowdose radiation depends on whether the cells have been stimulated to divide. The results show that significantly fewer breaks were observed in cells pretreated with 0.01 Gy in G 1 , S or G 2 than in those pretreated in G 0 . This suggests that the adaptation of cells to low doses requires mitogenic stimulation of lymphocytes. Concerning the lifespan of the adaptive response, data shows that lymphocytes treated with 1500 mGy at 40, 48, 66 h exhibited an adaptive response, and those treated later did not. Figure 6 shows the potential relationship between delta (δ) parameter and the time interval between priming and challenging dose (Δt). One can deduce that for approx. Δt > 100 h the Yonezawa effect completely disappears in human lymphocytes.
In the third analysed study (Shadley and Dai 1992), low doses were applied 12 h after culture (PHA application), at 18 h (first G 1 phase after PHA) after stimulation the higher dose was applied. In this case, however, DNA mutations (namely, the chromosomal aberrations) were analysed but the experimental investigation was quite similar to the previous two studies ) with a short T. Here, both aberration and deletion frequencies were reduced when a priming dose was applied, although there was a variability between the responses in samples from different donors. The results of this experiment show that G 1 lymphocytes are capable of exhibiting a cytogenetic adaptive response. Table 2 shows the exact raw data published in the three mentioned papers Shadley and Dai 1992) to estimate the values of α 0 , α 1 and α 2 parameters using both numerical methods. The two  . 6 The relationship between the delta (δ) parameter for DNA lesions and the time interval between priming and challenging dose (Δt) for human lymphocytes , where D 1 = 10 mGy, D 2 = 1.5 Gy, and T − Δt = 6 h. The dashed line represents the assumed potential linear trend (the best fit by δ = − 0.005 [h −1 ] Δt + 0.487 with R 2 ≈ 0.60) according to which the Yonezawa effect disappears above Δt ≈ 100 h. The straight line was selected as the simplest best fit according to the existing scatter of data points. Due to potential outliers, this fitting was also tested by the robust Bayesian regression method (Fornalski et al. 2010) but the result was practically the same (δ = − 0.005 [h −1 ] Δt + 0.483) first studies  can be simply treated as fully consistent and analysed jointly, because both of them show DNA lesions in the Yonezawa scheme. Thus, the results are: • for the SGA: 0 = 22.9 +0.5 −4.0 Gy −2 h −3 , 1 = 79.4 +5.5 −11.2 Gy −1 and 2 = 0.0832 +0.0093 −0.0082 h −1 . • for the L-BFGS-B algorithm: α 0 = 22.9 Gy −2 h −3 , α 1 = 79.5 Gy −1 and α 2 = 0.0832 h −1 .
The results obtained by both algorithms are qualitatively the same within the range of uncertainties, therefore both numerical methods work well and they can be treated on equal footing. In the case of the SGA, the uncertainties were estimated by the upper-lower bound method. Simulations were run for each worst case scenario assuming that numbers of lesions or mutations measured in experiments follow a Poisson distribution. As explained in the Appendix 2, in the L-BFGS-B algorithm it was not possible to calculate uncertainties of the parameters. Table 2 The source data used for estimation of α 0 , α 1 and α 2 parameters, selected for DNA lesions and chromosomal aberrations, based on three studies by Shadley and Collaborators (Shadley and Wolff 1987;Sahdley and Dai 1992)     In the next step, one can observe that the data and results of all three cited papers Shadley and Dai 1992) are fully comparable, regardless of chromosomal aberrations analysis in the last one. Therefore, they are able to be treated as a meta-data for joint analysis, which is carried out in the "Discussion".

Mutations
The analogical numerical research for mutations only can be investigated using Eq. (10). For this case, the comprehensive experimental data collected by Day et al. (2006Day et al. ( , 2007 were used (see Table 3) because they offer a large combination of situations to be used in numerical analysis. In the study by Day et al. (2007), Atm knockout heterozygous pKZ1 mice were treated by a whole-body X-irradiation with a priming dose (D 1 ) and, after a 4 h interval, a challenging dose (D 2 ). 3 days after exposure, the mice were sacrificed, their spleen and prostate removed and snap-frozen until needed for scoring of chromosome inversions. Both the 0.01 and 10 mGy priming doses caused a similar in magnitude adaptive response to the challenging dose of 1 Gy, although a single dose of 0.01 mGy seemed to induce a chromosomal inversion frequency by itself. Based on the results by Day et al. (2007), being a Atm knockout heterozygote (Hishiya et al. 2005) does not affect the response to single low radiation doses (used here) or the induction of an adaptive response for inversions.
The second analysed study (Day et al. 2006), used the very low priming doses of 0.001, 0.01, 1 or 10 mGy, which were administered 4 h before a challenge dose of 1 Gy. Mice were sacrificed 3 days after treatment to detect pKZ1 chromosomal inversion assay in the prostate. The results show that the 1-10 mGy priming doses reduced chromosomal aberration (inversion) frequency to below the spontaneous frequency level. The 0.001 mGy single priming dose had no significant effect, while the 0.01 mGy single priming dose increased the inversion frequency. Despite this, all four low and ultra-low priming doses caused an adaptive response when a 1 Gy challenging dose was administered (Day et al. 2006), which is presented in Table 3.

Discussion
The role of the time interval between two consecutive doses of ionizing radiation has been pretty well known for years and is the basis of e.g. dose fractionation. But the first experimental finding that this time interval is related to repair processes was conducted in 1940s by Prof. John R. Raper during his works in Manhattan Project (Raper 1947). At that time he irradiated groups of mice using beta radiation and observed a "recovery from the radiation damage" when the conditioning sublethal dose was applied prior to the test lethal dose. Raper's experiments were repeated a few years later and showed the same effect but no detailed explanation was proposed (Cronkite et al. 1950). Next, in the late 1950s, the historical Elkind-Sutton experiments were conducted (Elkind and Sutton 1959). Elkind and Sutton irradiated mammalian cells with two doses, which they called conditioning and test doses, analogically to Raper's terminology. They found that the time delay between doses can change Table 3 The source data used for estimation of α 0 , α 1 and α 2 parameters, selected for DNA mutations, based on studies by Day et al. (2006Day et al. ( , 2007  the shape of the survival curve of irradiated cells but both doses were large and no low conditioning dose was tested. Important changes appeared in 1980s and 1990s, when many studies began reporting that the radiation adaptive response effect showed strong results with a low priming dose and high challenging dose scheme (Olivieri et al. 1984;Shadley and Dai 1992;Yonezawa et al. 1990;Fan et al. 1990;Liu et al. 1992;Farooqi and Kesavan 1993;Cai et al. 1994). In the 1990s, this special case of adaptive response was called the priming dose effect or, a few years later, the Yonezawa effect (Wang et al. 2013(Wang et al. , 2018(Wang et al. , 2021Liu et al. 2019). Today, we can also call it the Raper-Yonezawa effect to account for some historical aspects surrounding it. Anyway, from the experimental point of view, the easiest way to test the adaptive response appearance is the priming dose scheme: the first small dose, called the priming or conditioning dose, and some time later the higher one, called challenging dose. That situation results in smaller detrimental effects than when the same big dose is given as a single dose -of course only in situations when the radiation adaptive response was activated.
In the XXI century there have been several biomathematical models which have tried to explain, among others, the priming dose effect (Schöllnberger et al. 2001;Yonezawa 2003, 2004;Esposito et al. 2010;Zhao and Ricci 2010;Wodarz et al. 2014;Devic et al. 2018Devic et al. , 2020Bondarenko et al. 2021), however, they are usually based on a set of differential equations, and in most cases they are not deeply related to the adaptive response phenomenon. The presented paper intends to explain the biophysical origin of the Yonezawa effect within the scope of a relatively simple model used so far in many of our papers (Fornalski et al. 2017;Dobrzyński et al. 2016Dobrzyński et al. , 2019. The considerations concern the effect of the use of a high challenging dose after a rather low priming dose given earlier. In the original papers by Yonezawa and coworkers (Yonezawa et al. 1990(Yonezawa et al. , 1996Matsubara and Yonezawa 2004;Matsumoto et al. 2004;Yonezawa 2003, 2004), originally for mice, the challenging dose was so high that it caused the systematic death of mice with time. However, when the priming dose was added, this process was substantially hindered. The proposed new model describing the Yonezawa effect contains only three free parameters. This model, however, is based on the radiation adaptive response probability function which makes it more grounded in biophysics. The model was validated on a group of experimental data from papers by Shadley and Dai 1992) and Day et al. (2006;. Those studies represent the most broad and comprehensive data in the Yonezawa scheme, which is necessary for proper parameter calculations to validate the model. Other exemplary papers, which presented the effect (Fan et al. 1990;Liu et al. 1992;Farooqi and Kesavan 1993;Cai et al. 1994), contain too few results to be valuable for the model's validation.
In practice, the model is reduced to a relatively simple equation, which connects the delta (δ) parameter ( Fig. 1) with priming and challenging doses, and the radiation adaptive response probability function. This equation, however, is presented in two versions: for DNA lesions, namely fast damage (Eq. 9), and mutations, which are rather later effects (Eq. 10). The quantitative results obtained by Eqs. (9) and (10) are very similar, therefore one can use one equation (Eq. 9) for both types of data to improve and unify the model's parameters, as well as to reduce some of uncertainties. We are aware that the number of lesions a long time after application of the challenging dose may only approximate the number of mutations. Therefore our model apparently simplifies the real situation in which the kinetics of late repairs should be better described.
We have to note that it is known (Berthel et al. 2019;Feinendegen et al. 2007;Mezentsev and Amudson 2011;Long et al. 2007;Ding et al. 2005) that low and high doses trigger separate groups of genes. As stated by Mezentsev and Amudson (2011), "Accumulating data suggest that the biological responses to high and low doses of radiation are qualitatively different", that may be linked to the activation of "radiation-responsive genes after high-and low-dose exposures". Therefore, our model in which the adaptive response at both doses is described by the same function of dose and time may not be adequate for a full explanation of the originally observed Yonezawa effect. This puts a natural limit on the value of the challenging dose. On the other hand, the alpha parameters including their uncertainties, which are connected with the biology of the studied object, endpoints, and proper timing of the whole process as well, may be so individual-dependent that the whimsical behaviour of Yonezawa effect can be expected.
The dose relationship presented by Eq. (9) has to be commented on here as well. Certainly, the right-hand side of Eq. (9) is not appropriate for very high doses D 2 because it tends to zero for that case. This problem is much wider: firstly, a very high value of D 2 would be lethal, thus no Yonezawa effect may be observed because of the organism's death. Secondly, as mentioned above, a very high value of D 2 activates different groups of genes responsible for DNA repair, which makes the radioadaptation more complicated and this is not reflected in our relatively simple model. Therefore, the presented approach has a very important limitation: the model works well when the challenging dose, D 2 , is not loo large.
The alpha parameters calculated for all the data in Table 2 are as follows: These results are related to human lymphocytes Shadley and Dai 1992) and their potential radioadaptation, which reaches its maximum for the priming dose of D 1 = 2/α 1 ≈ 25.2 mGy and Δt = 2/α 2 ≈ 24 h after its exposure. It is worth mentioning here that in the original Yonezawa and Smirnova mice model Yonezawa 2003, 2004), their theoretical consideration of the biological mechanisms responsible for the modifications of radiosensitivity is due to the change of radiosensitivity of the hematopoietic system. Consistently, the blood-forming system's cells were the objects of the studies (Smirnova and Yonezawa 2004).
The next finding, which can be observed in the detailed data from Shadley et al. studies (Shadley et al. 1987), is that the level of radioadaptation is connected with the cell cycle phase and its radiosensitivity. The analysis of the data shows that the strongest radiation adaptive response generated by the priming dose is connected with the lowest radiosensitivity of the cell, which seems to be quite natural (Fig. 7). It has long been known that cells in different cell cycle stages display different sensitivity to radiation. Cells in the late S phase are usually most radioresistant and cells in the M phase are most radiosensitive. A cellular response to DNA-damaging agents correlates not only with DNA replication and chromosome segregation stage. It also depends on the activation of the repair pathway maintaining the genetic integrity and activation of the cell cycle checkpoint. This complex mechanism is driven by the variety of key proteins concentrated in the cell (such as tumour suppressor proteins p53, inhibitor Gy −1 and 2 = 0.0845 +0.0085 −0.0060 h −1 . proteins p21, CHK kinases, ataxia telangiectasia-mutated ATM or ataxia-telangiectasia and RAD3-related ATR) and their inner cell signalling (Pawlik and Keyomarsi 2004). After a cell enters cycle arrest two major DNA doublestrand break repair pathways can occur: error prone NHEJ (Non-Homologous End Joining, dominant in G 1 /S phase) and HR (Homologous Recombination, present in late S or G 2 phase). The greater proportion of error-free repair in late S phase may explain its radioresistance (lower radiosensitivity). If the DNA damage is fully repaired, the cell cycle continues. The choice of which pathway becomes activated is determined by the conflict between maintenance and resection of the DNA ends (Maier et al. 2016).
The differences in radioadaptation are, however, relatively small but nonetheless, the adaptive response in phase G 0 -G 1 and after S (G 2 ) is smaller than in G 1 -S, which is clearly presented in Fig. 7. The mechanisms explaining such cell cycle effects on the adaptive response to ionizing radiation are not fully understood, but it has been suggested that the adaptive response phenomenon may be due to cell cycle changes (Cramers et al. 2005) or arrest (Syljuåsen 2019) caused by the low priming dose of radiation. Since sensitivity to radiation varies with cell cycle stage, changes in cell cycle distribution may be responsible for the radiation adaptive responses (Hafer et al. 2007). Recently, many studies have been focussing on pointing out which of the repair pathway components are taking majority in radioadaptation induction (Hafer et al. 2007;Hendrikse et al. 2000;Boothman et al. 1998).
The delta (δ) parameter, when the input data are too weak, can hardly be determined because less than 3 equations for 3 variables {α 0 , α 1 , α 2 } would give an infinite number of solutions. This is the reason why advanced numerical methods were necessary to calculate the model's parameters. The delta (δ) parameter, however, is very sensitive to even the smallest changes of alphas, which display the Yonezawa effect in some ranges only. This is clearly illustrated in Fig. 8, where the presented surface can easily reach maximal or minimal values of δ.
Indeed, the uncertainties of the alpha parameters inferred from the data shown in Table 2 indicate the high selectivity of the model used in the description of the Yonezawa effect. However, the calculated alpha parameters can be useful for the prediction of delta factor for other dose schemes, which can be planned in future experiments. Moreover, this method can be useful for the radiation adaptive response probability function assessment (like Fig. 7) for e.g. cell behaviour modelling during irradiation.
One has to note that the presented results were calculated for experimental cases where the Yonezawa effect was pretty well observed. But this is not always the case: this effect is observed in approx. 50% of all expected cases (Tapio and Jacob 2007), which means that the radiation adaptive response appearance is quite a selective process: the probability that this phenomenon Fig. 7 The non-normalized probability functions of radiation adaptive response in human lymphocytes  in phase G 0 -G 1 (blue solid line), in phase G 1 -S (orange dashed line) and after phase S (green dotted line), related to time (hours). The orange dashed line corresponds to the lowest radiosensitivity of the cell and therefore the strongest radioadaptation (color figure online) appears equals approx. 0.5, but when it does appear, it can be described by the p AR distribution with proper {α} parameters.
The last item to discuss is terminology: it is often found in scientific literature that the priming dose effect (which can be called the Yonezawa-or Raper-Yonezawa-effect) is practically equivalent to the radiation adaptive response phenomenon. However, this is not the case: the priming dose (Yonezawa) effect is just a special case of the adaptive response with a specific dose fractionation scheme. Another examples of the adaptive response is constant low-dose rate irradiation (Dobrzyński et al. 2016) or so called "radiation training" by many small dose pulses (Socol et al. 2020), see Appendix 1 ("Mutations in DNA"). Nevertheless, year by year we learn more about radioadaptation (UNSCEAR 2000;Tapio and Jacob 2007;Guéguen et al. 2019) therefore it is high time for its wide biophysical and mathematical description, at least of the most popular priming dose scenario.

Conclusions
The presented paper uses the radiation adaptive response theory to create a single equation (Eqs. (9) or (10)) for the Raper-Yonezawa (priming dose) effect (Fig. 1) when D 2 is not too high, and the adaptive response probability function after D 2 can be described by same formula as the one after the priming dose D 1 . The delivered equations were confronted with a group of experimental data on humans and mice to calculate the exemplary input parameters. For example 0 = 22.9 +0.5 −4.0 Gy −2 h −3 , 1 = 79.4 +5.5 −11.2 Gy −1 and 2 = 0.0832 +0.0093 −0.0082 h −1 are related to the human lymphocytes DNA lesions reduction in Yonezawa scheme. The relatively narrow range of these parameters indicates that their values may be strongly dependent on the individual case. Indeed, it is a set of parameters for one type of cells and their individual radiosensitivity. Therefore, the presented analysis shows that the level of radiation adaptive response is strictly connected with radiosensitivity (Fig. 7)-low radiosensitivity (so high radioresistance) is a result of a strong adaptive response of the cell, which is consistent with many recent findings presented in this paper. Finally, the proposed mechanistic model was quantified based on experimental datae.g. lesions in human lymphocytes and chromosomal inversions in mice-to be able to predict the Raper-Yonezawa effect for future experimental and theoretical investigations.

DNA lesions
Let us consider the relationship between repaired damage (lesions) and the adaptive response, dN = − N p AR dt, see Eqs. (6) and (2). The simplest solution of that equation, as an indefinite integral, is represented by Fig. 8 a The shape of delta (δ) parameter function from Eq. (9) for the exemplary sets of parameters and their ranges: α 0 = 36.21 Gy −2 h −3 , α 1 from 20 to 1200 Gy −1 , α 2 from 0.02 to 0.23 h −1 , D 1 = 10 mGy, D 2 = 1.5 Gy, Δt = 34 h; plot b represents the same situation but two parameters' ranges were narrowed: α 1 from 50 to 550 Gy −1 , and α 2 from 0.08 to 0.1 h −1 so one can calculate it as and assume such a solution as a function f (D,t). Thus, the definite integral has a form: For example, after some calculations and integration for N from N 0 to N(T) and for t from 0 to T, one can get This can be written in exact form as: where N 0 means the initial number of lesions (just after single D appearance, as described by Eq. (5)), N(T) is the actual remaining lesions, and T corresponds to the actual time. It should be noted, that Eq. (16) has an inverted Gompertzianlike shape in the function of time (Fornalski et al. 2020). Let us assume for further analyses, that the term is constant for a dedicated dose pulse, D, which is also a constant value in exact conditions. One can also note, that f(D,0) = ξ D . Additionally, for further considerations, let us calculate the probability connected with the appearance of the adaptive response from dose D after some time Δt: which is an analogical solution to Eq. (16).
After the next period of time, from Δt to T, the adaptive response from the same dose D is calculated as: However, when another dose is received in the moment Δt, it is necessary to use Eq. (16) with time T shifted by Δt: In the scenario of the typical Yonezawa effect, one can consider the three moments in time mentioned above: (1) the first one, let us call it a moment zero, t 0 = 0, when the first small priming dose (D 1 ) was received, (2) the second one after a short period of time, Δt > 0, when the challenging (much higher, D 2 > D 1 ) dose pulse was received, and (3) the third one, T, which corresponds to the actual moment, 0 < Δt < T, when the number of lesions is measured.
Next, one can note that in the moment of Δt, the number of damage (lesions) from dose D 1 which remain unrepaired equals to (see Eq. (15)): where N 0,1 denotes the initial number of damage (lesions) induced by dose D 1 in moment zero (assumed that N 0,1 = μ 0 + μ 1 D 1 , see Eq. (5)).
In time Δt, dose D 2 generates an additional number of damage, N 0,2 = μ 1 D 2 , where N 0,2 > N 0,1 (because D 2 > D 1 ). Therefore, in the moment of Δt, one needs to consider the total N 0,2 + N' 0,1 number of damage to repair. However, in the same moment, additional repair mechanisms from the second dose appear, which by assumption are additive to the first ones (because repair mechanisms from the first dose are still working). That way, in the later time, T, the number of unrepaired damage (lesions) equals: which can be denoted as N 1+2 for further considerations (see Eq. (7)). Let us note that as D 2 increases, the value of D2 also increases with D 2 2 , which causes a rapid decrease of N(T) (see analogical Fig. 3). To understand these relationships, Fig. 2 presents the time related probability of the adaptive response (p AR ) and the number of unrepaired Finally, at N 0,2 > 0, the main quantification of the Yonezawa effect, the parameter δ, in the moment of T equals: where N 2 = N 0,2 exp[f(D 2 ,T − Δt) − f(D 2 ,0)]. One can note, that the denominator in the left-hand-side term of Eq. (23) represents the situation where the single dose D 2 (without a priming dose) is given in the time t 0 + Δt, as shown in Fig. 2.
The parameter δ can vary from 1 to some minimal value δ min which can even be lower than zero (but in most studies it is assumed that δ min = 0). Generally, to keep the delta parameter above δ ≥ 0, one needs to fulfill the condition of which guarantees that the adaptive response signal from dose D 1 will be still working when challenging D 2 appears. The Yonezawa effect completely disappears when the time distance (Δt) between D 1 and D 2 is too large for the activation of repair mechanisms of the priming dose. Parameter δ attains its minimal (negative) value of which is a clear information that in such a case the Yonezawa effect disappears.

Two priming doses
The two priming doses case is presented in Fig. 4. After the first time interval (i.e. just before the moment of Δt 1, so just before D 2 appears) one can use reasoning analogous to the previous case: Eq. (18) can be rewritten as − ∫ Δt 1 0 p AR dt , which can be used analogically to the Eq. (22). After the second time interval (i.e. just before the moment of Δt 1 + Δt. 2, so just before D 3 appears), Eq. (19) (29)). The process of mutations repair can be written as: The parameter δ can be now defined as 1 -N 1+2+3 /N 3 , so: which describes the Yonezawa effect when the challenging dose D 3 is applied after two priming doses D 1 and D 2 . Equation (33) is equivalent to the Eq. (11).

Multiple priming doses
To present more general situation, let us consider n identical priming doses D separated by the same time shift Δt. The situation where several single dose pulses are applied to an organism is called a "radiation training" (Socol et al. 2020). The challenging dose D* (where D* > D) is applied Δt after the last of n doses D were applied. Using the same reasoning as above, one can present the general equation for the parameter δ as: One can consider the situation where the number of dose pulses is infinite (n ⟶ ∞) and the time distance between them tends to zero (Δt ⟶ 0). This situation is equivalent to chronic irradiation by a constant dose rate (Ḋ = const), which can be easily found in areas with high natural background radiation (Dobrzyński et al. 2015a(Dobrzyński et al. , 2015b, for example. However, Eq. (34) would give an infinite level of adaptive protection ( lim Δt→0;n→∞ = 1 ), which is not possible.
To avoid that problem, one should note that in this particular situation the probability function of the radiation adaptive response (Eq. (3)) saturates at some constant value: which was originally calculated in (Dobrzyński et al. 2016), where Ḋ = const. This result will modify Eq. (16) to N(T) = N 0 exp(− P C ) = const which allows one to make analogical calculations as in previous subchapters. This is, however, not the subject of this paper, because chronic irradiation with constant dose-rate is not considered in the Yonezawa effect. Still, one can pose a question whether the natural radiation can function as a factor reducing possible effects of eventual higher doses (Mortazavi et al. 2012). (32b) To conclude, the radiation adaptive response effect (namely, radioadaptation) is a wide phenomenon. Special case of radioadaptation is the priming dose effect (called the Raper-Yonezawa effect). Another example of the adaptive response is e.g. constant low-dose rate irradiation (see Eq. (35)), like in high background radiation areas (Dobrzyński et al. 2015a(Dobrzyński et al. , 2015b.

Simplified Genetic Algorithm (SGA)
The proposed methodology starts with the well-known least squares method, where Eqs. (9) (or (10)) has defined the function to be fitted. This function is complicated, thus an iterative method of finding the minimum of sum of squared residuals had to be used. Thus, the most suitable method of parameters {α} assessment is the genetic algorithm. However, the use of the classical genetic algorithm (Banzhaf et al. 1998;Whitley 1994) is not suitable for finding a global minimum, for several reasons. One of them is limited possibility of coding a broad range of real numbers (such as parameters {α}) into artificial "genes" with keeping sufficient resolution. Another is that two steps of the classical genetic algorithm (namely crossover and mutation) have a tendency to change these parameters rapidly, resulting in an escape from the best solution neighbourhood. Significant changes to the algorithm were introduced including removal of the crossover step and applying small random changes directly to the parameters instead of randomly mutating artificial "genes". A much simpler yet more effective algorithm mixing the features of genetic algorithms with particle swarm optimization 2,3 and simulated annealing 4 was then obtained (Fig. 9).
Simulation starts with creating a set ("a population") of 99 individuals. Each individual represents all three-parameter values: α 0 , α 1 and α 2 , chosen randomly as real positive numbers at the beginning of simulation. These values are used to calculate the fitness function of each individual, which in our case is reciprocal of the sum of squared residuals (each residual is the difference between δ measured in the experiment and δ calculated with given values of α 0 , α 1 and α 2 ). Next, a new set of individuals is selected from the old one with the probability of being chosen proportional to the fitness function of a given individual. This way individuals with the best fitness can be duplicated and these with the worst can be eliminated. Next, small random changes are applied to parameters in the new population, which becomes the current one. These steps make up one loop ("a generation") of the algorithm, which is repeated until an end condition is fulfilled. The condition can be defined as being unable to find better parameters as the simulation progresses or simply reaching the maximum generation number, which in our case was 6,000,000.
Since there is a risk of not finding a global solution in a single simulation, simulations conducted by the described algorithm were run at least 50 times and manually checked for obvious mistakes (i.e. when the algorithm got stuck in a local solution or wasn't able to even get on the right path to global solution). These mistakes can be spotted by comparing the fitness functions of many simulations. In our case, it was safe to assume that if they differ more than a 1% from the best fitness function obtained in any simulation, the result should be ignored. Luckily, there were not many mistakes (less than 10%). Moreover, correct results were used to define the region of interest and more precise values of parameters were obtained by a brute-force search. 5

Bounded Limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS-B) algorithm
The second approach also uses χ 2 function to find the best estimation of alpha parameters. In this case, however, the optimization procedure uses a quasi-Newtonian algorithm called BFGS (Broyden-Fletcher-Goldfarb-Shanno) using a limited amount of computer memory with additional bound constraints of variables (Byrd et al. 1995;Zhu et al. 1997). L-BFGS-B method uses an estimate of the inverse Hessian matrix to drive its search in the variable space. The advantages of this algorithm are fast convergence and relatively low computational complexity. However, one should also be aware of its shortcomings: the L-BFGS-B algorithm may turn out to be divergent when the starting point is far from the solution sought and it does not guarantee finding the global minimum-it may happen that the parameters found only correspond to the local minimum. To minimize the risk that the found parameters would correspond to a local minimum instead of a global minimum, the algorithm was started from multiple points in the parameter space within a reasonable range. The L-BFGS-B algorithm should not be used for very flat functions. The limitations of numerical precision make it impossible to compute the gradient of such a function. This creates the last disadvantage: the calculated uncertainties of the fitted function's parameters are inconsistent, thus the uncertainties were calculated by SGA method only.