Quantitative modeling of multigenerational effects of chronic ionizing radiation using targeted and nontargeted effects

Stress response signals can propagate between cells damaged by targeted effects (TE) of ionizing radiation (e.g. energy depositions and ionizations in the nucleus) and undamaged “bystander” cells, sometimes over long distances. Their consequences, called non-targeted effects (NTE), can substantially contribute to radiation-induced damage (e.g. cell death, genomic instability, carcinogenesis), particularly at low doses/dose rates (e.g. space exploration, some occupational and accidental exposures). In addition to controlled laboratory experiments, analysis of observational data on wild animal and plant populations from areas contaminated by radionuclides can enhance our understanding of radiation responses because such data span wide ranges of dose rates applied over many generations. Here we used a mechanistically-motivated mathematical model of TE and NTE to analyze published embryonic mortality data for plants (Arabidopsis thaliana) and rodents (Clethrionomys glareolus) from the Chernobyl nuclear power plant accident region. Although these species differed strongly in intrinsic radiosensitivities and post-accident radiation exposure magnitudes, model-based analysis suggested that NTE rather than TE dominated the responses of both organisms to protracted low-dose-rate irradiation. TE were predicted to become dominant only above the highest dose rates in the data. These results support the concept of NTE involvement in radiation-induced health risks from chronic radiation exposures.

www.nature.com/scientificreports/ valuable because they produce data that span wide ranges of dose rates and radiation types (e.g. mixtures of gamma, beta and alpha-emitting radionuclides) applied over many generations. Analysis of such data, along with those from laboratory experiments, can enhance our understanding of ionizing radiation effects, including NTE.
Here we used a mechanistically-motivated mathematical model of TE and NTE, based on our previous work 41,44,45,48,49 , to analyze published embryonic mortality data for plants (Arabidopsis, Arabidopsis thaliana) 50 and rodents (bank vole, Myodes or Clethrionomys glareolus) from the Chernobyl nuclear power plant accident region 51 . These data are valuable because they were collected relatively soon after the disaster (starting within 1-2 years after the accident in 1986), when radiation dose rates were highest and their effects were likely to be most prominent, and covered many generations of studied plants and animals, with each subsequent generation being exposed to a lower dose rate due to radionuclide decay and ecological processes. Our goal in analyzing such data was to investigate the model-predicted roles of NTE versus TE as function of radiation dose rate and time since the start of exposure, and to compare these results across species. We believe that such studies provide potentially useful evidence about the mechanisms of radiation-induced health risks at low doses/dose rates.

Materials and methods
Mathematical model. We previously developed a simple mathematical model of radiation-induced NTE and implemented it on several human and animal data sets 41,44,45,48,49 . Although the signaling pathways involved in NTE are complex and incompletely understood, their consequences can be quantitatively modeled by using the following set of assumptions: (1) irradiated cells "activate" other cells in an "on-off " (binary) manner by NTE signals, causing them to enter into a prolonged stressed state, which can be transmitted across generations (e.g. epigenetically). (2) "Activated" cells accumulate damage at an elevated rate. (3) Eventually they can revert to the background state, but this process may be very slow. (4) NTE-induced damage adds to the damage produced by direct traversal of targets by radiation (TE).
The adaptation of this modeling approach to the data sets on embryonic mortality during multi-generational exposure to chronic radiation, where the dose rate decreases over time, is described by the following system of ordinary differential equations. Here R is the radiation dose rate, t is time after the start of irradiation, P a is the average probability of cells to be in an "activated" state due to NTE signals induced by radiation, and Y is radiation-induced yield of the damage of interest (embryonic mortality in this case).
The parameters k 1 (dose −1 ) and c 3 (time −1 ) represent cell "activation" and deactivation by NTE, respectively. This notation is the same as in our original papers on NTE modeling 48,49 . The parameter k bac (time −1 ) represents background Y formation, k TE (dose −1 ) represents Y formation by TE, k NTE (time −1 ) represents the proportionality constant relating P a to Y, and κ (time −1 ) represents Y removal (e.g. due to damage repair and selection against heavily damaged cells/organisms in the population).
Without radiation exposure (R(t) = 0), the system (Eqs. 1, 2) approaches an equilibrium condition where P a (t) = 0 and Y(t) = k bac /κ. For irradiation at constant dose rate (R(t) = R c ), the system can be solved analytically as follows: Assuming that all model parameters are greater than zero, after long times at constant dose rate the system tends towards the following equilibrium solution: However, here we are interested in the situation where radiation dose rate is not constant, but decreases exponentially over time since the Chernobyl nuclear accident according to the following equation, where R 0 is the initial dose rate right after the accident and λ is its rate of decrease: www.nature.com/scientificreports/ Of course, this exponential time dependence for the dose rate is only a rough approximation for environmental exposures after an accident because complex processes of radionuclide migration and bioaccumulation in different soil layers and biomass of different organisms are not explicitly accounted for. However, the exponential dependence does account for the main trend in dose rate after an accident with the minimal number of adjustable parameters.
When the exponential dose rate dependence on time (Eq. 9) was substituted into the system of model equations (Eqs. 1, 2), we could no longer find an analytic solution and solved the system numerically using the dsolve numeric procedure in Maple 2019 software. Since Y(t) in our model represents the mean number of lethal events (a continuous number > 0), whereas the observed data are probabilities of embryonic mortality (bounded between 0 and 1), we used the following formula based on an assumed Poisson distribution of lethal events to generate model predictions for the probability of embryonic mortality (P mort ), which is the probability of observing ≥ 1 lethal event: This approach was used to fit predicted P mort values to the plant and rodent data sets, as described below.
Data sets. The plant (Arabidopsis thaliana) data set was obtained from Tables 1-8 in the publication by Abramov et al. 50 . Three variables were used for analysis: (1) time since the beginning of irradiation (since the Chernobyl accident in this case), (2) radiation dose rate at each time point at each studied location, and (3) the corresponding embryonic mortality yield at each time point at each studied location.
The radiation dose rate at each of 14 locations (Chernobyl, Shepelichi, Stechanka, Tolstyi Les, Yanov-0, Yanov-1, Yanov-2, Yanov-3, Yanov-4, Yanov-6, Yanov-7, Yanov-8, Yanov-10, Damba) was calculated based on the "monthly dose" for the sum of γ and β radiations reported by Abramov et al. 50 , converted to units of µGy/h (assuming 1 month = 30 × 24 h). Dose rate was ln-transformed to bring the error distribution closer to Normal. Robust linear regression with time since 1986 (in years) as the independent variable, and ln-transformed dose rate (in µGy/h) as the dependent variable, was performed (using the rlm function in R 4.0.2 software) separately for each location. This regression generated location-specific estimates of the dose rate at time zero (parameter R 0 ) and the exponential rate of decrease for the dose rate (λ, years −1 ). For those locations where there were not enough time measurements to perform the regression (e.g. dose rate measured only at one time), λ was set to approximately zero (10 -3 years −1 , to avoid a singularity in solving the model equations) and R 0 was set to the mean of available dose rate measurements. The robust procedure was selected instead of ordinary least squares regression to minimize the potential effect of "outlier" data points on the R 0 and λ estimates. The embryonic mortality yield was based on the "Chimeras for lethals, %" values reported by Abramov et al. 50 , with the addition of "the maximum frequency observed in control plants (5%)". The resulting data set composed of 35 data points is provided in Table 1.
The rodent (Clethrionomys glareolus) data set was obtained from Table 3 in the publication by Ryabokon et al. 51 . The radiation dose rate at each of 3 locations (called sites 2, 3 and 4) was calculated based on the "Wholebody absorbed dose rate" for the sum of γ and β radiations reported by Ryabokon et al. 51 , converted from µGy/ day to µGy/h. For two instances of background conditions where no dose rate value was reported by Ryabokon et al. 51 , we assumed 0.05 µGy/day = 0.00208 µGy/h.
As for the plant data described above, dose rate was ln-transformed to bring the error distribution closer to Normal. Robust linear regression was performed separately for each location to estimate parameters R 0 and λ. For those locations where there were not enough time measurements to perform the regression, λ was set to 10 -3 years −1 and R 0 was set to the mean of available dose rate measurements, for the reasons described above. The embryonic mortality yield was based on the "Mean embryonic lethality" values (totals for before and after implantation) reported by Ryabokon et al. 51 . The resulting data set composed of 12 data points is provided in Table 2.
Both Arabidopsis thaliana and Clethrionomys glareolus have short generation times, generally less than 1 year, and sometimes 2 per year, under the studied conditions 50,51 . Therefore, these data sets analyzed here extend to multiple generations of each species.
Model fitting procedure. The model equation for embryonic mortality yield under an exponentially decreasing dose rate (Eq. 9 substituted into Eq. 2) was fitted separately to each data set (plant and rodent) using a customized iterative algorithm implemented in maple 2019 software. This algorithm is described in Supplementary Methods online. Uncertainties of the best-fit model parameters on each data set were estimated by generating 1000 randomly perturbed versions of each data set by bootstrapping with replacement. The fitting procedure was applied to each perturbed data set, starting with the best-fit parameter combination previously found on the original (not perturbed) version of the corresponding data set. The 2.5th and 97.5th percentiles of the distribution of best-fit values of each parameter over the 1000 perturbed data sets were used as estimates of the 95% confidence interval (CI) of the selected parameter.

Results
Our model-based analysis was able to describe the main patterns of the analyzed plant and rodent data sets using a small number of adjustable parameters ( These results suggest that NTE rather than TE dominated the responses of both organisms to protracted low-dose-rate irradiation (Figs. 3, 5). In both species, the TE contribution to the radiation response because numerically substantial, relative to NTE, only at dose rates ≥ 10,000 µGy/h. Such high dose rates were present in (10) www.nature.com/scientificreports/ the plant data set (Table 1), but not in the rodent data set ( Table 2). For rodent data, the TE component therefore represents an extrapolation. However, because the dose rate range ≥ 10,000 µGy/h was investigated by many laboratory studies, the importance of TE for intense radiation exposures is expected. The presumed importance of NTE at lower dose rates in the current analysis is based on the properties of the model used here (described in "Materials and methods" section). The model attributes nonlinear concave radiation response shapes to NTE, and the observed nonlinearities in the analyzed data at low dose rates are fitted by NTE parameters. Model parameters for both NTE and TE, and their uncertainties, are shown in Table 3 for plants and rodents.
Since the model consists of a system on nonlinear differential equations, parameters can affect the fit in combinations (e.g. k 1 /c 3 , k bac /κ), and the influences of each parameter or combination are not the same. For these reasons, some parameters had much larger uncertainties than others (Table 3). For example, parameter k 1 , which represents cell activation by non-targeted effects, was very uncertain in both data sets.
Nevertheless, comparison of best-fit radiation response parameter values for plants and rodents resulted in the following observations: (1) Parameter c 3 , which represents cell deactivation from the NTE state, and parameter κ, which represents radiation-induced damage removal over time, were much smaller for rodents, than for plants.
(2) Parameter k NTE , which represents the relationship between NTE activation probability and damage formation, was similar across data sets. (3) The TE parameter (k TE ) for rodent data is probably unreliable because high dose rates at which TE are likely to become important were not represented in the data set. www.nature.com/scientificreports/ These differences in parameter values (particularly c 3 and κ) resulted in qualitative differences in predicted   www.nature.com/scientificreports/ radiation responses for plants are rodents. Specifically, in situations of exponentially decreasing dose rate, radiation-induced embryonic mortality was predicted to continue increasing for considerable time even though dose rate was monotonically decreasing. In Arabidopsis, the predicted time since the start of irradiation needed for embryonic mortality to peak was a few years/generations, and a rapid decline was predicted afterwards (Fig. 4).   www.nature.com/scientificreports/ In voles, the predicted time needed for the effect to peak was much longer than the maximum time in the data set (10 years after the accident, corresponding to ~ 20 generations) (Fig. 6). In summary, model-based analysis suggested that both Arabidopsis and voles are susceptible to presumably NTE-mediated transgenerational effects such as embryonic mortality, but in voles the NTE state is much more persistent and the damage is removed from the population much more slowly. Consequently, in Arabidopsis the peak yield of radiation-induced damage was predicted to occur within a few years/generations after the start of exposure to an exponentially decreasing dose rate (Fig. 3), whereas in voles the peak damage yield could occur several decades (> 20 generations) after the start of exposure (Fig. 6), when the dose rate was already very low.

Discussion
Since individual cells in a multicellular organism constantly interact by various forms of signaling, damage caused by ionizing radiation in some cells generates signals that affect the responses of other cells and can induce them to enter into a prolonged stressed state. This state can be transferred even across generations (probably epigenetically) 11,34 . This phenomenon implies that descendants of severely irradiated individuals can continue to have elevated rates of damage, even if they themselves were exposed much less severely or not at all. The data sets analyzed here 50,51 represent likely examples of this phenomenon under environmental-not laboratory-conditions. In particular, the vole data show elevation of embryonic mortality over multiple generations, even though radiation dose rates during this period were decreasing 51 . This pattern is strongly suggestive of non-targeted effects (NTE) such as genomic instability 51 . Our simple mathematical model was able to generate behaviors consistent with this finding by being capable of generating predictions where the radiation effect continues to grow even when the dose rate is reduced dramatically from its initial maximal value.
Radiation-induced NTE were previously detected and studied in Arabidopsis thaliana 18,19,22 and in rodents [52][53][54] , including transgenerational increases in tumorigenesis 55 and chromosome aberrations 34 . These findings provide indirect support for the plausibility of the conclusions of the current modeling analysis of Arabidopsis and vole embryonic mortality data from the Chernobyl accident region 50,51 , which suggest a very important role for NTE in responses to protracted low dose rate irradiation in both species. Although Arabidopsis and rodents are phylogenetically very distant and differ strongly in genome size and intrinsic radiosensitivities, both appear to be susceptible to radiation-induced NTE. Transgenerational damage, attributed mainly to NTE by our radiation response model (described in "Materials and methods" section), was particularly persistent in voles, but was removed from the population more rapidly in Arabidopsis.
In voles, this phenomenon resulted in continuously increasing rather than decreasing embryonic mortality over 10 years (~ 20 generations) after the accident, even though radiation dose rates decreased dramatically over this period 51 . In Arabidopsis, the response was maximal a few years/generations after the accident, and decreased afterwards 50 . These qualitative differences in terms of presumed long-term persistence of NTE could be caused by a variety of factors, such as differences in the mode of reproduction between animals (direct production of gametes) and vascular plants (alternation of gametophyte and sporophyte generations).
Both types of radiation response patterns-transiently peaking or persistent transgenerational effects-were decently described by our simple TE + NTE model, with different parameters for the plant and rodent data sets. The NTE-dominant interpretation of the radiation response behaviors at low dose rates is consistent with other studies of these and similar data 46,51,56 . However, alternative explanations are also possible, e.g. changes in radiation damage repair efficiency as function of dose rate 50 . The modeling approach and data sets also had important limitations in terms of model simplicity, fitting methodology, and limited data sample size and quality. For example, it is not known which specific NTE mechanisms could be involved in the observed radiation effects on Chernobyl animals and plants. Many environmental variables, such as temperature, humidity, soil type and resource availability, which probably had big effects on the studied organisms, were not available in the analyzed data sets and could not be modeled explicitly, resulting in relatively low R 2 values for model fits. Despite these limitations, we believe that the results obtained here support the concept of NTE involvement radiation-induced health risks from chronic exposures.