Mechanistic Modeling of Dose and Dose Rate Dependences of Radiation-Induced DNA Double Strand Break Rejoining Kinetics in Saccharomyces cerevisiae

Mechanistic modeling of DNA double strand break (DSB) rejoining is important for quantifying and medically exploiting radiation-induced cytotoxicity (e.g. in cancer radiotherapy). Most radiation-induced DSBs are quickly-rejoinable and are rejoined within the first 1–2 hours after irradiation. Others are slowly-rejoinable (persist for several hours), and yet others are essentially unrejoinable (persist for >24 hours). The dependences of DSB rejoining kinetics on radiation dose and dose rate remain incompletely understood. We hypothesize that the fraction of slowly-rejoinable and/or unrejoinable DSBs increases with increasing dose/dose rate. This radiation-dependent (RD) model was implemented using differential equations for three DSB classes: quickly-rejoinable, slowly-rejoinable and unrejoinable. Radiation converts quickly-rejoinable to slowly-rejoinable, and slowly-rejoinable to unrejoinable DSBs. We used large published data sets on DSB rejoining in yeast exposed to sparsely-ionizing (electrons and γ-rays, single or split-doses, high or low dose rates) and densely-ionizing (α-particles) radiation to compare the performances of the proposed RD formalism and the established two-lesion kinetic (TLK) model. These yeast DSB rejoining data were measured within the radiation dose range relevant for clonogenic cell survival, whereas in mammalian cells DSB rejoining is usually measured only at supra-lethal doses for technical reasons. The RD model described both sparsely-ionizing and densely-ionizing radiation data much better than the TLK model: by 217 and 14 sample-size-adjusted Akaike information criterion units, respectively. This occurred because: the RD (but not the TLK) model reproduced the observed upwardly-curving dose responses for slowly-rejoinable/unrejoinable DSBs at long times after irradiation; the RD model adequately described DSB yields at both high and low dose rates using one parameter set, whereas the TLK model overestimated low dose rate data. These results support the hypothesis that DSB rejoining is progressively impeded at increasing radiation doses/dose rates.

Introduction statistically robust comparisons of kinetic models. DSB rejoining has been measured by several methods: e.g. neutral filter elution, pulsed-field gel electrophoresis, or by surrogate markers such as γ H2AX foci [22,29,[42][43][44]. In mammalian cells, the first two methods produce reliable results only at supra-lethal radiation doses (generally > 20 Gy), at which cells remain metabolically functional for some time, but are clonogenically dead [22,29,31]. The third method is applicable to lower doses, but the kinetics of foci accumulation and decay can be quite different from the underlying DSB rejoining kinetics. In Saccharomyces cerevisiae, which has a much smaller genome than mammalian cells, DSB rejoining can be quantified directly at doses relevant for clonogenic cell survival [45].
Extensive data sets on DSB rejoining in S. cerevisiae were produced by Frankenberg-Schwager et al. [38][39][40][41]. Petite mutant yeast (diploid strain 211 Ã B) lack mitochondrial DNA, which facilitates radio-labeling of nuclear DNA. The yeast cells were lysed on top of a neutral sucrose gradient (5-20%) and the released DNA was sedimented by centrifugation [38]. Gradients were fractionated onto glass fiber filters, dried, and treated with a toluene-based scintillation liquid. The percentage of total radioactivity of DNA in each fraction as a function of the sedimented volume yielded DNA profiles which were used to determine the number of radiationinduced DSBs. The average number of DSBs per molecular mass of DNA was calculated by computer simulation of random breakage as applied to the DNA of unirradiated cells and by fitting these calculated curves to the DNA profiles obtained from irradiated cells [38].
This methodology was applied by Frankenberg-Schwager et al. to study DSB rejoining after sparsely-ionizing or densely-ionizing radiation under non-growth conditions, which allowed rejoining to occur, but prevented confounding of the results by cell proliferation. Here, we analyzed data sets for the following irradiation scenarios: 1. Sparsely-ionizing 30 MeV electrons, delivered at a high dose rate (7800 Gy/h), either as a single dose (300-2400 Gy) followed by rejoining time of 0-72 h, or as split doses (900-2400 Gy/dose) separated by a 16-48 h interval and followed by rejoining time of 0-24 h [38]. These data were taken from Figures 3-7 of reference [38]. The linear energy transfer (LET) was approximately 0.21 keV/μm.
For convenience, we classified the data as high dose rate (HDR, data sets 1 and 3) or low dose rate (LDR, data set 2) exposures. The data points, which were presented graphically in the cited publications by Frankenberg-Schwager et al. were digitized using GetData Graph Digitizer software (http://www.getdata-graph-digitizer.com/). In some cases (particularly for data set 3), not all experimentally measured data points were published: instead, only the mean and standard error were reported. In this situation we used a conservative assumption that the number of data points was three, and assigned values to these points so that the reported mean and standard error would be reproduced: one point was set equal to the mean, and the other two were symmetrically positioned around it at an appropriate distance to mimic the standard error. The total number of data points in the three data sets analyzed here was 278: 186 for sparsely-ionizing and 92 for densely-ionizing radiations. The values of these data points are provided in the S1 Appendix.

Mathematical models
In the proposed RD model (schematically summarized in Fig 1), radiation (with dose rate R) produces three DSB classes (DSB 1 , DSB 2 , DSB 3 ) with yields (per unit dose) of k 1 , k 2 and k 3 , respectively. Quickly-rejoinable DSBs (DSB 1 ) and slowly-rejoinable DSBs (DSB 2 ) are rejoined with rates v 1 and v 2 , respectively; DSB 3 are unrejoinable. Radiation converts DSB 1 to DSB 2 and DSB 2 to DSB 3 with rates proportional to parameters q 1 and q 2 , respectively. The model contains seven parameters: k 1 , k 2 , k 3 , v 1 , v 2 , q 1 and q 2 . These assumptions are mathematically In the RD model, radiation (with dose rate R) produces three DSB classes (DSB 1 , DSB 2 , DSB 3 ) with yields (per unit dose) of k 1 , k 2 and k 3 , respectively. DSB 1 and DSB 2 are rejoined with rates v 1 and v 2 , respectively; DSB 3 are unrejoinable. Radiation converts DSB 1 to DSB 2 and DSB 2 to DSB 3 with rates proportional to parameters q 1 and q 2 , respectively. In the TLK model, radiation also produces three DSB classes with yields (per unit dose) of c 1 , c 2 and c 3 , respectively; DSB 1 and DSB 2 are rejoined with rates λ 1 and λ 2 , respectively; DSB 3 are unrejoinable.DSB 1 and DSB 2 are "fixed" to become DSB 3 with rates ε 1 and ε 2 , respectively. DSB 1 interact with each other (by quadratic mis-rejoining) with rate η 1 , and DSB 2 interact with each other with rate η 2 . DSB 1 interact with DSB 2 with rate η 1,2 . Details are described in the main text. represented by the following system of differential equations: These equations are analytically solvable, and the solutions are provided in the S2 Appendix for each DSB class at time T after a single radiation dose D is delivered with dose rate R. Linear quadratic (LQ) approximations to these solutions (which are computationally convenient and can be compared to analogous approximations of other kinetic models [46,47]) are also included in the S2 Appendix.
In the TLK model [34] (schematically summarized in Fig 1), radiation also produces three DSB classes with yields (per unit dose) of c 1 , c 2 and c 3 , respectively; DSB 1 and DSB 2 are rejoined with rates λ 1 and λ 2 , respectively; DSB 3 are unrejoinable. DSB 1 and DSB 2 are "fixed" to become DSB 3 with rates ε 1 and ε 2 , respectively. DSB 1 interact with each other (by quadratic mis-rejoining) with rate η 1 , and DSB 2 interact with each other with rate η 2 . DSB 1 interact with DSB 2 with rate η 1,2 . The model contains ten parameters: c 1 , c 2 , c 3 , λ 1 , λ 2 , ε 1 , ε 2 , η 1 , η 2 , and η 1,2 . These assumptions are represented by the following system of differential equations: Here, we have included a term for direct yield of unrejoinable DSBs (parameter c 3 ) for completeness, because although this term was absent from the original TLK model [34], it was introduced in subsequent similar formalisms [48]. There is no analytical solution to Eq 2.

Model fitting procedure
The RD and TLK models were fitted to data by maximizing the log likelihood, using optimization routines in Maple 171 software. All sparsely-ionizing radiation data sets (HDR, single and split doses, and LDR, single doses) were fitted together to assess how well can each model describe all of these types of exposures using one set of parameters. Although two distinct types of sparsely-ionizing radiation (30 MeV electrons and 60 Co γ-rays) were used to produce these data sets [38,39], their biological effects were likely to be similar. Densely-ionizing radiation data (HDR, single doses only) [40,41] were fitted separately. All parameters were restricted to ! 0 to maintain mechanistic plausibility.
We assumed a Gaussian error distribution with constant variance for all data points. This assumption was reasonable because the errors introduced during measurement of the data (DSB/cell yields) were not estimated explicitly [38] and were likely to be the same for all data points. The log likelihood function (LL M,i ) under assumption of constant variance, is described by the following equation: Here, the index M represents either the RD or the TLK model, i represents the data set (either sparsely-ionizing or densely-ionizing radiation), N (i) is the total number of data points for the i-th data set, P M,i,(j) are predictions of the M-th model at the j-th data point, and O i,(j) are observed data values. The constant term ln(2π)+1 was included for completeness, but had no effect on the comparison of model performances and on parameter estimation.
Analytic solutions for the RD model were used to calculate P M,i,(j) , and the likelihood function (Eq 3) was optimized by a sequential quadratic programming (SQP) algorithm [49] implemented in Maple 171 software. The probability of finding the global maximum (rather than local maxima) was enhanced by using 100 random initial conditions for the model parameters. For the TLK model, only numerical solutions could be obtained. Consequently, a customized optimization procedure (described in the S2 Appendix), also with 100 random initial conditions, written in Maple 171 was used for fitting the TLK model to data.
Absolute goodness of fit (GOF) for the analyzed models under assumption of constant variance was assessed by exploratory data fitting using three methods: (1) visual inspection of model fits and the data; (2) calculation of the coefficient of determination, R 2 ; (3) linear regression of model predictions P M,i,(j) vs. data points O i,(j) . Model fits were assumed to have no gross systematic deviations from the data when the 95% confidence intervals (CIs) for the regression of model predictions vs. data included 0 for the intercept and 1 for the slope.
Additional exploratory calculations showed that an alternative assumption of error magnitudes proportional to O i,(j) values, which resulted in replacement in Eq 3 of the term , reduced GOF for both the RD and TLK models because the data points at low doses, which had small values of O i,(j) , affected the fit more strongly than other points.

Estimation of model parameter uncertainties
Uncertainties (95% CIs) for best-fit model parameter values were estimated by profile likelihood [50] as follows: 10,000 Monte-Carlo-generated parameter values in the vicinity of the best-fit values were used to estimate the critical contour of the log likelihood function, which is based on the asymptotic X 2 behavior of the log likelihood distribution.

Information theoretic model selection
Ranking of models by relative support from the data, taking into account sample size and number of parameters, can be based on the Akaike information criterion with sample size correction (AICc). AICc has gained popularity for this purpose in various fields [51,52]. For comparing non-linear models, AICc is preferable to methods that rely on reduced X 2 or R 2 [53][54][55][56]. One of the most useful features of AICc is that it allows the evidence for structurally distinct models to be compared, without the need for models to be "nested" or to belong to the same class. The model that loses the least amount of Kullback-Leibler information relative to other compared models has the lowest AICc value and is considered as best-supported by the data at hand. The equation for AICc for the M-th model on the i-th data set (AICc M,i ) is given below, where K M is the number of adjustable parameters and LL M,i is the maximized log-likelihood value (calculated using Eq 3): The likelihood of the M-th model relative to other tested model(s), called the evidence ratio (ER M,i ), can be expressed as: Here, AICc min,i is the lowest AICc value generated by the set of models being compared. If ΔAICc M,i > 6, then the evidence ratio ER M,i becomes < 0.05, suggesting that the M-th model has much poorer support from the data than the best-supported model.
In contrast to the AICc, the Bayesian information criterion BIC [57] does not adjust for sample size effects [52]. The F-test requires arbitrary assumptions about test type (forward, backward or stepwise) and about the test significance threshold (α-value) [52], and provides no readily-interpretable information on relative support for each model from the data. Consequently, we believe that AICc is the most useful and convenient tool for selecting among plausible mechanistic models using DSB rejoining data.

Comparison of model performance
Best-fit predictions of the RD and TLK models after single-dose exposures are compared with the data in Figs 2-5. Visual inspection clearly shows that at long times (! 24 h) after sparselyionizing or densely-ionizing radiation only the RD model reproduced the upwardly-curving dose responses for remaining unrejoined DSBs, whereas the TLK model predicted linear dose responses and underestimated the data. For example, the mean measured number of DSBs/cell 72 hours after 2400 Gy of HDR 30 MeV electrons was 30.2 (range: 24.7-35.6). The corresponding best-fit prediction from the RD model was 28.4, whereas the TLK model predicted 19.5 (Fig 2). The mean measured number of DSBs/cell 72 hours after 600 Gy of α-particles was 20.1 (range: 9.6-30.9). The corresponding best-fit prediction from the RD model was 17.0, whereas the TLK model predicted 14.5 ( Fig 5). Therefore, the tendency of the TLK model to underestimate the data at long times and high doses was evident for both sparsely-and densely-ionizing radiations, although it was more clear with the former. The number of split-dose data points (10) was too small for robust conclusions (i.e. both models visually approximated these data reasonably), and consequently we did not show the fits to these data graphically.
Although the number of HDR data points (157) for sparsely-ionizing radiation was much larger than the number of LDR points (29), the RD model reasonably described DSB yields at both high and low dose rates using one set of parameters (Fig 4). In contrast, the TLK model was unable to do so: its fit was dominated by HDR data and LDR data were strongly overestimated (Fig 4). For example, the mean measured number of DSBs/cell 22 hours after 1650 Gy of LDR γ-rays was 5.6 (range: 4.2-6.4). The corresponding best-fit prediction from the RD model was 4.0, whereas the TLK model predicted 16.8 (Fig 4).
These visually-apparent differences in model performance were quantified by relative and absolute GOF assessments (Tables 1 and 2). The RD model described the sparsely-ionizing radiation data dramatically better than the TLK model: the sum of squared deviations from the data was 3.1-fold lower for the RD model, and this translated into a massive difference in support of 217 AICc units (Table 1).
On the smaller α-particle data set, which included only HDR exposures, the difference in model performances was reduced, but the RD model was still strongly favored: by 14 AICc units (Table 1). R 2 was higher and 95% CIs for the intercept and slope of linear regression of model predictions vs. observed data point values were often narrower for the RD model vs. the TLK model for HDR data ( Table 2). Systematic deviations of predictions from LDR data were suggested for the TLK model, for which the regression slope CIs did not include unity, but no such deviations were detected for the RD model ( Table 2).
The AICc scale being logarithmic, the relative support from the sparsely-ionizing data for the TLK model vs. the RD model (i.e. the evidence ratio, Eq 5) is exp[-217/2], which is zero for all practical purposes. For densely-ionizing radiation data, the TLK model evidence ratio is exp [-14/2], or approximately 1/1400. These marked differences in support are dominated by the differences in the log likelihoods (Eq 3) of the RD and TLK models (105 units for sparsely-ionizing and 3.5 units for densely-ionizing radiation data), rather than by a mere penalization for extra parameters. In fact, even if the TLK model were artificially favored, for example, by fixing three of its parameters (η 1 , η 2 , and η 1,2 ) at best-fit values, thereby reducing the number of adjustable parameters to the same number as in the RD model (7), the RD model would remain the best-ranked model by a still very comfortable margin of 210 AICc units for sparsely-ionizing and 7.1 units for densely-ionizing radiation data.
However, it is important to note that the TLK model performed much better than its predecessor-the LPL formalism [37]. When applied to DSB rejoining data, the LPL model contains two classes of DSBs (rejoinable and unrejoinable) and 4 adjustable parameters: the yields of each DSB class per unit dose, and coefficients for linear and quadratic DSB rejoining processes. On sparsely-ionizing radiation data, this simpler model had dramatically less support than the TLK and RD models (by 314 and 531 AICc units, respectively). R 2 for the LPL model on HDR singledose data was only 0.25, whereas it was 0.91 and 0.95 for the TLK and RD models, respectively ( Table 2). On densely-ionizing data, the LPL model was worse than the TLK and RD models by 204 and 218 AICc units, respectively, and its R 2 was 0.52 (vs. 0.96 for the TLK and RD models).
These results suggest that the TLK formalism, which contains three DSB classes, represents an important improvement over the LPL model with two classes. However, the RD model, which contains radiation-induced conversion of DSB classes, yields a further major improvement for the yeast data at hand. The improved performance of the RD model relative to the TLK model is caused specifically by the dose/dose rate dependence of conversion of DSB classes: if the conversion is allowed to occur without influence by radiation (i.e. the terms q 1 ×R and q 2 ×R in Eq 1 are replaced with q 1 and q 2 , respectively), model performance plummets by 210 and 7 AICc units for sparsely-ionizing and densely-ionizing radiation, respectively. In other words, without radiation-induced conversion of DSB classes, the RD model performs equivalently to the TLK model.

RD model parameter values
The best-fit RD model parameter values (Table 1) suggested that sparsely-ionizing radiation directly induced approximately 8 times more quickly-rejoinable than slowly-rejoinable DSBs (based on the parameter ratio k 1 /k 2 ). The direct yield of unrejoinable DSBs was not detectable: parameter k 3 had a best-fit value of zero. Consequently, the suggested interpretation is that unrejoinable DSBs were produced only indirectly by conversion of slowly-rejoinable DSBs with a rate proportional to the radiation dose rate and to parameter q 2 . Conversion of quicklyrejoinable DSBs (with a rate proportional to the radiation dose rate and to parameter q 1 ) was an important indirect source of slowly-rejoinable DSBs. The rejoining rates for these two DSB classes differed by a factor of 24 (based on the parameter ratio v 1 /v 2 ).
In contrast, the RD model suggested that densely-ionizing α-particles directly produced only slowly-rejoinable DSBs; whereas the direct yields of quickly-rejoinable and unrejoinable DSBs were not detectable: parameters k 1 and k 3 had best-fit values of zero ( Table 2). The yield of slowly-rejoinable DSBs per unit dose was approximately 18.5-fold higher than for sparselyionizing radiation (based on the ratio of parameter k 2 values). However, their rejoining rate (parameter v 2 ) was the same for both sparsely-ionizing and densely-ionizing radiation ( Table 1). The rate (per unit dose) of conversion of slowly-rejoinable DSBs to unrejoinable ones was approximately 2-fold higher for α-particles vs. sparsely-ionizing radiation (Table 1).

Discussion
We analyzed large published data sets on DSB rejoining in yeast (Saccharomyces cerevisiae) exposed to sparsely-ionizing and densely-ionizing radiation [38][39][40][41] to enhance the understanding of dose/dose rate dependences of radiation-induced DSB rejoining kinetics. Specifically, we tested the hypothesis that the fraction of slowly-rejoinable/unrejoinable DSBs increases with increasing dose/dose rate. This hypothesis was implemented in a new RD model, whose performance was compared to that of the established TLK model [34].
The strengths of the current study include rigorous information theoretic comparison of the performances of two models with different assumptions about the mechanisms of DSB Table 2. GOF assessment for the RD and TLK models. R 2 = coefficient of determination. The intercept (i) and slope (s) of linear regression of model predictions vs. observed data point values were used to assess GOF. Systematic deviations of predictions from the data (marked in bold font) are suggested if the 95% CIs (shown in parentheses) of the intercept do not include zero and/or if the 95% CIs of the slope do not include unity.  Table 1. Information theoretic performance assessment (by ΔAICc) and best-fit parameter values for the RD and TLK models. ΔAICc = 0 suggests the strongest support from the data for the given model among all tested models, whereas ΔAICc > 6 suggests poor support. Parameter meanings are provided in the caption to Fig 1 and in the main text.

Radiation Model ΔAICc
Best-fit parameter values (95% CIs) rejoining, using extensive data sets which include different radiation types, doses, dose rates, and rejoining times. The weaknesses include the use of S. cerevisiae (rather than mammalian cell) data. Yeast data are useful because DSB rejoining was measured in the dose range relevant for clonogenic cell survival [45], but some aspects of DSB rejoining in yeast and mammalian cells are different, e.g. because S. cerevisiae relies on homologous recombination to a greater extent than many other organisms [58,59], and chromatin structure is also different in yeast and mammals [60]. The numbers of radiation-induced DSBs per cell at the mean lethal dose for yeast and mammalian cells are comparable [60,61]. However, because the yeast genome is much smaller than mammalian genomes [62], radiation doses used for yeast studies (! 300 Gy of sparsely-ionizing radiation) are much higher than those relevant for clonogenic survival of mammalian cells (< 20 Gy). This difference in absolute dose magnitudes may have an effect on DSB rejoining kinetics which is independent of the number of DSBs per cell, but dependent on other factors such as the yields of reactive oxygen and nitrogen species (ROS and RNS). For example, radiation-induced oxidants reduce chromosomal movements by affecting the meiotic actin cytoskeleton, thereby potentially preventing unfavorable chromosome interactions [63].

Sparsely
In addition, at radiation doses relevant for yeast survival, cellular enzyme inactivation and changes in cell membrane integrity and permeability may play more important roles in cell death [64], than at the much lower doses relevant for mammalian cell survival. Despite these potential drawbacks, we believe that the current study provides potentially clinically-relevant insight into mechanistic modeling of DSB rejoining. In particular, because the TLK and RD models are based on different mechanistic assumptions, comparison of their performances on the same data enabled us to make the following observations about which putative mechanisms are most (or least) useful for describing DSB rejoining kinetics.
The TLK model contains 3 parameters (η 1 , η 2 , and η 1,2 ) for quadratic interactions between different DSBs within/between DSB classes (Fig 1). If some/all of these parameters have positive values, the model structure (Eq 2) implies that, as the dose of HDR radiation (and hence the yield of DSBs just after irradiation) increase, the rate of DSB removal increases because quadratic interactions between DSBs become more frequent. Consequently, the quadratic interactions mechanism predicts that the fraction of DSBs which remain unrejoined at a given time after irradiation should decrease with increasing radiation dose. The dose response for unrejoined DSBs at a given time after irradiation is, therefore, predicted to be downwardlycurving (i.e. to have a negative second derivative).
In contrast, the data analyzed here suggest the opposite pattern: at higher doses, the fraction of DSBs which remain unrejoined at each measured time after irradiation increases, rather than decreases (Figs 2, 3 and 5) [61]. Consequently, the observed dose response for unrejoined DSBs at a given time after irradiation is upwardly-curving (i.e. has a positive second derivative). This disparity between predictions and observations explains why the best-fit values of parameters η 1 , η 2 , and η 1,2 in the TLK model became zero (Table 1): when quadratic interactions between DSBs do not occur, the TLK model predicts a linear dose response for unrejoined DSBs (Figs 2 and 5), which is closer to the data than a downwardly-curving one.
The mechanism of DSB fixation, which is also present in the TLK model, could not be detected in the data analyzed here: the best-fit values of parameters ε 1 and ε 2 were zero ( Table 1). If either/both of these parameters were restricted to positive values, the fit quality did not change substantively, and the main outcome was reduction of the best-fit value of direct induction of unrejoinable DSBs (parameter c 3 ). In other words, the TLK model relied on unrejoinable DSBs (produced directly by radiation and/or by fixation of rejoinable DSBs) to describe the HDR data for both sparsely-ionizing and densely-ionizing radiation (Table 1). If this mechanism was turned off by setting parameter c 3 to zero, TLK model performance strongly worsened: by 7.6 AICc units for sparsely-ionizing and by 17.7 units for densely-ionizing radiation.
An important drawback of the TLK model's reliance on unrejoinable DSBs to approximate HDR data was that the best-fit yield of such DSBs (parameter c 3 ) was too high to describe LDR data, leading to overestimation of these data (Fig 4). In other words, in the context of the TLK model, the observed DSB data at long times after HDR sparsely-ionizing radiation exposure appeared to be inconsistent with the data after LDR exposure.
These results are in agreement with previous analyses which suggested that the RMR model does not describe yeast DSB rejoining data well [65]. The current study shows that the TLK model, which is more flexible than the RMR model because it allows for different DSB classes, achieved decent absolute GOF (Table 2). However, the TLK model performed very poorly relative to the RD model on the data at hand because it lacked the proposed mechanism of radiation-induced conversion of DSB classes.
In contrast, the proposed mechanism of radiation-dependent conversion of DSB classes present in the RD model (Fig 1) allowed the upwardly-curving dose response shape for DSBs left unrejoined after HDR radiation to be reproduced (Figs 2 and 5). There was also no inconsistency between HDR and LDR data in the context of the RD model (Figs 2-4): both data sets were described with adequate GOF (Table 2) using one set of best-fit parameters. These results were obtained thanks to the fact that the RD model predicts that, at increasing doses/dose rates, DSB rejoining is increasingly impeded because more and more quickly-rejoinable DSBs are converted to slowly-rejoinable DSBs, some of which are in turn converted to unrejoinable ones.
The qualitative differences in dose and dose rate dependences of DSB rejoining predicted by the RD and TLK models (using best-fit parameter values) can be visualized in Fig 6. The TLK model predicts that the yield of unrejoinable DSBs per unit dose is independent of both dose and dose rate, whereas the yield of slowly-rejoinable DSBs depends on dose rate but not on dose. The RD model, however, predicts that the yields of both slowly-rejoinable and unrejoinable DSBs depend on both dose and dose rate (Fig 6). Our analysis suggests that the proposed mechanism of radiation-dependent conversion of DSB classes strongly enhanced the ability of kinetic models to describe DSB rejoining in yeast. In contrast, the analyzed data did not support the notion (represented in the TLK model) that interactions between different DSBs constitute an important route of DSB removal. However, interactions between DSBs which lead to transition into a different class of DSB would still be compatible with the data. Our analysis also suggests that the time needed to rejoin a particular DSB depends mainly on the properties of this DSB and/or on radiation damage to nearby chromatin. These factors can be represented by the generic term "DSB complexity", which depends not only on radiation type, but also on dose and dose rate.
The previously mentioned GLOBLE model [17,28] provides a conceptually similar interpretation and also predicts that the fraction of slowly-rejoinable/unrejoinable DSBs should increase with radiation dose. Specifically, the GLOBLE model assumes that the genome is organized into chromatin loops, and that the presence of one or multiple DSBs inside a given loop is associated with different damage severity and rejoining rates. The GLOBLE approach has been applied to mammalian cell data, but it is not clear whether its assumptions are applicable to yeast because chromatin organization in yeast differs from that in mammals [66,67].
In contrast, the RD model can be viewed as a more generalized approach where the nature of DSB classes and the mechanisms of conversion between them are not specifically tied to chromatin structure. For example, when repair enzymes encounter a complex DSB, and/or are oxidatively damaged during the repair process (which is more likely to occur at higher radiation doses/dose rates), rejoining of the DSB may become stalled. Conceptually similar explanations were considered by the authors of the data sets analyzed here, but were never implemented as a mathematical model [68]. Such explanations, as well as those based on chromatin structure, are consistent with the RD formalism.
The LQ approximation to RD model solutions for the number of DSBs which remain unrejoined at long (effectively infinite) rejoining times (Eqs. A3-4 of the S2 Appendix) provides additional insight into dose rate effects within the context of the RD model. At all dose rates, the linear DSB dose response component is composed of the yield of unrejoinable DSBs (parameter k 3 ). The quadratic component is composed of the product of the yield of slowly-rejoinable DSBs (parameter k 2 ) and the conversion coefficient for slowly-rejoinable to unrejoinable DSBs (parameter q 2 ). This quadratic component is absent at low dose rates. Thus, lowering the dose rate reduces the quadratic dose response term, but the mechanistic interpretation of this phenomenon is not the same as in other kinetic models: in other models the quadratic term is reduced because, at low dose rates, most DSBs are rejoined before they can interact with each other; whereas in the RD model the quadratic term is reduced because, at low dose rates, most quicklyrejoinable DSBs are rejoined before radiation can convert them to slowly-rejoinable ones.
These conclusions are important for mechanistic understanding and quantitative modeling of DSB rejoining kinetics and have potential clinical relevance, e.g. for optimizing cancer radiotherapy. However, because we cannot exclude the possibility that the results of this study are specific to S. cerevisiae, additional testing of the RD model needs to be performed on mammalian cell data. Specifically, it remains to be determined what molecular mechanisms define the differences between DSB classes and explain radiation-induced conversion between them.