Modeling the Effect of Hypoxia and DNA Repair Inhibition on Cell Survival after Photon Irradiation

Mechanistic approaches to modeling the effects of ionizing radiation on cells are on the rise, promising a better understanding of predictions and higher flexibility concerning conditions to be accounted for. In this work we modified and extended a previously published mechanistic model of cell survival after photon irradiation under hypoxia to account for radiosensitization caused by deficiency or inhibition of DNA damage repair enzymes. The model is shown to be capable of describing the survival data of cells with DNA damage repair deficiency, both under norm- and hypoxia. We find that our parameterization of radiosensitization is invariant under change of oxygen status, indicating that the relevant parameters for both mechanisms can be obtained independently and introduced freely to the model to predict their combined effect.


Introduction
Radiation therapy is one of the cornerstones of cancer care, where~50% of patents receive radiation during the course of disease [1]. Radiobiological modeling is an integral part of radiation oncology and radiation therapy, used to predict normal and tumor tissue response, which is of particular significance when moving towards personalized radiation treatment, including treatment gap corrections, normal tissue tolerance predictions, optimization of therapy determined by predictive assays, multi-modality schedule design, and the simulation of clinical trials [2]. In this work, we present the first step towards development of our modeling platform called "UNIfied and VERSatile Engine" (UNIVERSE), within which we aim to integrate multiple biological responses and mechanisms relevant for describing radiation response of different cell types. In this manuscript, we focus on describing the cellular response of a particularly radioresistant tumor sub-population-hypoxic cells. Cells under hypoxic conditions exhibit increased radioresistance, and tumors containing hypoxic regions have significantly worse chances of successful treatment with radiotherapy [3][4][5]. The increased radiosensitivity of cells in the presence of free oxygen is usually explained by the oxygen fixation hypothesis: Molecular oxygen has the ability to react with radicals produced in the DNA, thereby fixating the damage, preventing the direct chemical restoration of the DNA radical by reacting with H + [3]. The concept of the oxygen enhancement ratio (OER) is classically used to quantify the dependence of cell survival on the oxygenation status. The OER is usually defined as the ratio between the doses needed to induce the same survival fraction in a hypoxic and a normoxic environment [6]. One of the key strategies to overcome the radioresistance in hypoxic cells is a dual treatment, i.e., a combined treatment with photon irradiation and administration of radiosensitizing drugs, such as DNA damage response (DDR) inhibitors [7][8][9][10].
To model the response of hypoxic tumor cells to dual treatment, we based our approach on a previously published model of our group [11]. Similarly to models described by other groups [12][13][14], the here presented model divides the cells' nucleus into equally sized subvolumes containing about 2 Mbp, dubbed giant loops [15][16][17]. Those giant loops containing either exactly one or two and more double strand breaks (DSB) are classified as isolated DSB (iDSB) or complex DSB (cDSB), respectively. Giant loops have been identified as possible critical targets, inside which multiple lesions resist swift repair [18][19][20][21]. It was shown that computed numbers of iDSB (N iDSB ) and cDSB (N cDSB ) matched well with observed frequencies of quickly (iDSB) and slowly (cDSB) repaired DSB in rejoining studies [22,23]. The core model of our choice [11,13] associates both classes of lesions with so-called lethality parameters, K iDSB and K cDSB , respectively. These parameters are the corresponding probabilities for lesions of a given class to become lethal, meaning the cell loses its potential to further proliferate. Complex lesions are considered to be significantly more lethal than isolated lesions, as each of them poses a high risk for chromatin loss [14] and they remain unrepaired for a prolonged time [24][25][26]. It is indeed found that K cDSB is several magnitudes larger than K iDSB in cases derived from experimental data [11,26].
Regarding the hypoxic cell population, Carlson et al. [6] made compelling arguments for an interpretation of the OER as the ratio of doses needed to induce the same total amount of DSBs in hypoxic and normoxic cells. They further suggested the replacement of the OER term with the name hypoxia reduction factor (HRF O 2 DSB ) at a given oxygen concentration [O 2 ] within this context [27]. In line with the above ideas, we introduced the HRF O 2 DSB as a parameter into our model, which solely modifies the initial total yield of DSB (N tDSB ). Following evidence that the oxygen status has no effect on the DSB rejoining rates [6], the lethality parameters were assumed to be invariant under the change of oxygenation. We could show in our previous work, that the sole introduction of the HRF O 2 DSB parameter was sufficient to describe survival data from literature and that the derived HRF O 2 DSB values were well described by a parameterization suggested by Carlson et al. [6] (a function of oxygen concentration [O 2 ], inspired by the initial studies of Alper and Howard-Flanders [28]) [11]. In the work presented here, we investigate whether our model is capable of describing an experimental set of cell survival data containing five cell lines at three different oxygen levels each, and whether the derived HRF O 2 DSB values are in accordance with the formerly introduced parameterization.
To further extend the model and describe the particular case of hypoxic tumor cells response to dual treatment, radiotherapy, and DDR inhibition, we introduced a so-called radiosensitization factor (RSF) that modifies K iDSB . This is based on observations by Hufnagl et al. [26], that the increased radiosensitivity of repair-deficient cell lines could be accounted for by increasing the lethality parameter of isolated DSB, while keeping the lethality parameter for complex DSB constant. Complex DSB are argued to pose such severe challenge to the DDR that any change in the repair capabilities of the cell has no effect on their lethality parameter. To validate this extension, we benchmarked our model in two scenarios, using experimental data obtained from cells in which one of the two key radiation-induced DDR molecules, DNA-dependent protein kinase (DNA-PK) or ataxia-telangiectasia mutated (ATM), was impaired [29]. First, we studied the robustness of the model to predict survival data of DNA-PK-deficient mutants of CHO cells. Secondly, we tested the model to predict the survival of two human lung cancer cell lines with pharmacologically-inhibited ATM. Both scenarios were investigated under normoxia and hypoxia.

Modeling Hypoxia-Induced Radioresistence
Our survival data for five different cell lines (A549, H460, H1437, B16, Renca) exposed to three distinct oxygenation levels (normoxia 20% [O 2 ], 1% [O 2 ], 0.5% [O 2 ]) were fitted using our model ( Figure 1). K iDSB and K cDSB were fitted for each cell line as free parameters to the normoxic data. Keeping these parameters fixed, for each cell line, the specific HRF O 2 DSB for the two hypoxic conditions were fitted. The determined numerical values of each parameter can be found in Table 1. Our model shows excellent capability in describing the acquired data, indicating consistency with our earlier work [11]. Furthermore, the HRF O 2 DSB values derived from our data were in accordance with the parametrization of HRF O 2 DSB as a function of oxygen concentration published in [11], as shown in the bottom right panel of Figure 1.

Modeling Hypoxia-Induced Radioresistence
Our survival data for five different cell lines (A549, H460, H1437, B16, Renca) exposed to three distinct oxygenation levels (normoxia 20% [ ], 1% [ ], 0.5% [ ]) were fitted using our model ( Figure 1). and were fitted for each cell line as free parameters to the normoxic data. Keeping these parameters fixed, for each cell line, the specific for the two hypoxic conditions were fitted. The determined numerical values of each parameter can be found in Table 1. Our model shows excellent capability in describing the acquired data, indicating consistency with our earlier work [11]. Furthermore, the values derived from our data were in accordance with the parametrization of as a function of oxygen concentration published in [11], as shown in the bottom right panel of Figure 1.

Modeling Cell Survival of DNA-PK-Impaired Cell Lines
Survival data of the CHO cell line and two of its DNA-PK response-deficient mutants (V3 cell line is DNA-PKcs-deficient and xrs-5 cell line is Ku80-deficient) under normoxia and hypoxia were gathered from Cartwright et al. (Figure 2) [30]. The cell line-specific parameters K iDSB and K cDSB were found by fitting our model to the normoxic survival data of the wild-type cells. The HRF O 2 DSB of the cell line was then derived by fitting the model to the hypoxic survival data, keeping the aforementioned lethality parameters (K iDSB and K cDSB ) fixed. By increasing only the lethality parameter of isolated DSB (K iDSB ) by a radiosensitization factor (RSF) for each of the mutant cell lines, their survival under normoxia was able to be fitted accurately. The numerical values of the parameters can be found in Table 2. By applying the HRF O 2 DSB derived from the wild-type data and the RSF derived for each mutant under normoxia, the survival of the two mutant cell lines under hypoxia could be predicted satisfactorily. However, deviations can be observed at the highest reported doses. Survival data of the CHO cell line and two of its DNA-PK response-deficient mutants (V3 cell line is DNA-PKcs-deficient and xrs-5 cell line is Ku80-deficient) under normoxia and hypoxia were gathered from Cartwright et al. (Figure 2) [30]. The cell line-specific parameters and were found by fitting our model to the normoxic survival data of the wild-type cells. The of the cell line was then derived by fitting the model to the hypoxic survival data, keeping the aforementioned lethality parameters ( and ) fixed. By increasing only the lethality parameter of isolated DSB ( ) by a radiosensitization factor (RSF) for each of the mutant cell lines, their survival under normoxia was able to be fitted accurately. The numerical values of the parameters can be found in Table 2. By applying the derived from the wild-type data and the RSF derived for each mutant under normoxia, the survival of the two mutant cell lines under hypoxia could be predicted satisfactorily. However, deviations can be observed at the highest reported doses.

Modeling Cell Survival of Cell Lines with Pharmacologically-Inhibited ATM
The same approach was applied to our data containing two of the initially presented cell lines (H460 and H1437) exposed to different concentrations of an ATM inhibitor (ATMi) and irradiated under normoxia and hypoxia ( Figure 3). First, the cell line-specific lethality parameters ( and ) were derived by fitting our model to the data of cells irradiated under normoxia and without drug treatment. Second, the of each cell line was derived by fitting the model to the data of cells receiving no drug but irradiated under hypoxia, keeping the lethality parameters and fixed. Third, RSFs for the lethality parameters of the isolated lesions ( ) were derived for each cell line and drug concentration by fitting the model to the data of cells at each drug concentration under normoxia. The numerical values of the derived parameters can be found in Table  3. Again, by applying the derived from the non-treated cells and the RSF values found for each drug concentration irradiated under normoxia, the survival of the two cell lines exposed to the combination of different drug concentrations and hypoxia could be predicted very well.

Modeling Cell Survival of Cell Lines with Pharmacologically-Inhibited ATM
The same approach was applied to our data containing two of the initially presented cell lines (H460 and H1437) exposed to different concentrations of an ATM inhibitor (ATMi) and irradiated under normoxia and hypoxia ( Figure 3). First, the cell line-specific lethality parameters (K iDSB and K cDSB ) were derived by fitting our model to the data of cells irradiated under normoxia and without drug treatment. Second, the HRF O 2 DSB of each cell line was derived by fitting the model to the data of cells receiving no drug but irradiated under hypoxia, keeping the lethality parameters K iDSB and K cDSB fixed. Third, RSFs for the lethality parameters of the isolated lesions (K iDSB ) were derived for each cell line and drug concentration by fitting the model to the data of cells at each drug concentration under normoxia. The numerical values of the derived parameters can be found in Table 3. Again, by applying the HRF O 2 DSB derived from the non-treated cells and the RSF values found for each drug concentration irradiated under normoxia, the survival of the two cell lines exposed to the combination of different drug concentrations and hypoxia could be predicted very well.

Discussion
Our model provided an excellent description of the survival data of five cell lines and three oxygen levels presented in Figure 1, confirming the applicability of the previously published framework [11]. It underlines the hypothesis introduced in our former work, that cell survival under hypoxic conditions can be described in a first approximation by keeping the defined lethality of isolated and complex lesions invariant and only modifying the overall induction of DSBs by a given factor ( ). In contrast to our former publication, we determined the lethality parameters and by fitting them to the experimental data. This practice might lead to better predictions, as up to now, these values were recalculated from provided or fitted LQ model parameters ( and values), based on a Taylor expansion at low doses of our model equations. However such approximate recalculations remain to be crucial in cases in which only the and values are available but not the full set of cell survival data. Furthermore, we were able to show that the derived values from our data coincided well with the parametrization as a function of oxygen concentration introduced in our former publication. The highest value for both oxygen levels was obtained from A549 cells. Including these data points, the mean of the derived values deviated by 0.08 and 0.19 for 1% [ and 0.5% [ , respectively. If one excludes the A549 data, the mean of the derived values only deviated from the prediction by 0.01 and 0.07 for 1% [ and 0.5% [ , respectively. This is further evidence for this parameterization to be a widely applicable estimate for the in cases where the data to derive the exact value are not available.
The idea of increasing the lethality of isolated lesions under the presented model in order to describe the increased cell killing observed for repair deficient cell lines has already been expressed and successfully demonstrated by Hufnagl et al. within the GLOBLE model [26]. They consider the

Discussion
Our model provided an excellent description of the survival data of five cell lines and three oxygen levels presented in Figure 1, confirming the applicability of the previously published framework [11]. It underlines the hypothesis introduced in our former work, that cell survival under hypoxic conditions can be described in a first approximation by keeping the defined lethality of isolated and complex lesions invariant and only modifying the overall induction of DSBs by a given factor (HRF O 2 DSB ). In contrast to our former publication, we determined the lethality parameters K iDSB and K cDSB by fitting them to the experimental data. This practice might lead to better predictions, as up to now, these values were recalculated from provided or fitted LQ model parameters (α and β values), based on a Taylor expansion at low doses of our model equations. However such approximate recalculations remain to be crucial in cases in which only the α and β values are available but not the full set of cell survival data. Furthermore, we were able to show that the derived HRF O 2 DSB values from our data coincided well with the HRF O 2 DSB parametrization as a function of oxygen concentration introduced in our former publication. The idea of increasing the lethality of isolated lesions K iDSB under the presented model in order to describe the increased cell killing observed for repair deficient cell lines has already been expressed and successfully demonstrated by Hufnagl et al. within the GLOBLE model [26]. They consider the lethality of complex lesions K cDSB as being fixed, as each complex lesion poses a "significant burden for the cell", irrespective of the DNA damage repair capabilities of the cell. Further, they argue, that in NHEJ deficient cell lines, the lethality of isolated lesions increases depending on the cell cycle status of the cells. While we adapt and fully agree with the notion of a fixed lethality of complex damages independent of the repair capabilities of a cell, we had no information on the cell cycle distributions underlying the data taken from Cartwright et al. [30] and found different radiosensitizing factors (RSF) for the two DNA-PK response-deficient mutants analyzed. Therefore we decided to introduce the RSF as a free parameter in UNIVERSE, fitted to each mutant cell line. The RSF values found for the two CHO mutants V3 and xrs5 indicate that the probability of an isolated lesion to become a lethal lesion increases through DNA-PK response deficiency in these cell lines by a factor of about 10 and 15, respectively. Even though both RSF factors are related to the functional DNA-PK repair activity, the difference in the RSF of both cell lines might be retraced to the fact, that both are deficient of different enzymes taking part in the DNA-PK response: While xrs-5 cells are deficient of the Ku80 DNA-PK subunit [30,31], V3 cells lack the catalytic subunit for DNA-PK [30,32]. It is, however, unclear if and how the RSF is coupled to the activity of both proteins. Furthermore, different extents of remaining expression or compensation by other repair proteins might also lead to the observed difference in the RSF values of both mutants. One would have to compare several groups of cell lines with the exact same deficiencies to gain a better understanding of the underlying dependencies. More importantly, we could show that the survival of the mutant cell lines under hypoxia were well predicted by our model by combining the HRF O 2 DSB determined based on the wild-type data and keeping the RSF values derived from the normoxic data of each mutant invariant.
Since both DNA-PK and ATM are essential molecules for irradiation-induced DNA damage repair, recruited to the DNA damage sites [33][34][35], we assumed that pharmacological inhibition of ATM should also lead to an increased lethality of isolated damage sites, similarly to the genetic models in which DNA-PK is deficient, with a concentration-dependent effect. Indeed, the observed survival in cells exposed to different concentrations of the ATM inhibitor under normoxia was described with high fidelity by introducing an RSF for isolated lesions. Further, the values of the derived RSF illustrate the increasing lethality of isolated lesions with increasing drug concentrations, leading to stronger inhibition of the repair, while staying below the value derived for a pathway deficiency. It is also to be noticed that the RSF values for both H460 and H1437 cells are fairly similar. However, more cell lines have to be analyzed to investigate whether this can be extended to a general trend and if the presented method is generally applicable to other repair-inhibiting drugs. Nevertheless, we could show that the RSF values derived from normoxic data can accurately describe the survival in hypoxic conditions by introducing the HRF O 2 DSB derived from cells without drug treatment. Taken together, we could show that impairment of DNA damage repair, both in repair-deficient cell lines and cells treated with a DDR inhibitor, could be accounted for by a manipulation of the lethality of isolated lesions with an RSF in our model. Moreover, this modified lethality could stay invariant under change of oxygen supply, while sustaining good predictive capabilities. Thus, the uniqueness of our approach lies in its capability to describe two separate cell response mechanisms in any combination using minimal input parameters, which can be separately derived. Based on this, we believe that our approach has a high potential to implement further cellular mechanisms in order to produce predictions tailored to diverse clinical scenarios.

Experimental Data from Literature
Experimental data used to benchmark the model on survival of genetically DDR-deficient cell lines were taken from [30]. The data used to benchmark the model on survival of NCI-H460 (H460) and NCI-H1437 (H1437) cells in which the DDR was pharmacologically inhibited using an ATM inhibitor under normoxia and hypoxia, as well as the survival data of A549 cells under hypoxia and normoxia, were taken from [10].

Cell Culture, Clonogenic Survival Assay, and Irradiation
For the validation of our model under hypoxic conditions, additional experiments were performed using Renca (murine renal carcinoma; American Type Culture Collection, Manassas, VA, USA) and B16-Blue ISG (murine melanoma; Invitrogen, Thermo Fischer, Waltham, MA, USA) (B16) cells. Both cell lines were grown in RPMI 1640 Medium (Gibco, Thermo Fischer, Waltham, MA, USA) supplemented with 10% fetal bovine serum (FBS)(Merck Millipore, Darmstadt, Germany) at 37 • C and 5% CO 2 atmosphere. Experiments in hypoxic conditions were performed at 0.5 or 1% O 2 and 5% CO 2 using a custom hypoxic chamber (C-chamber; Biospherix, Parish, NY, USA), including an online monitoring controller for O 2 and CO 2 concentrations (ProOx and ProCO2 model; Biospherix, Parish, NY, USA). Fifty cells per well in a 96-well format were seeded not more than 16 h before irradiation. Hypoxic irradiation was performed after incubation for 4 h under respective oxygen conditions. Cells were irradiated in the sealed hypoxia chamber with a dose series of photons of 1, 2, 4, or 8 Gy and thereafter incubated under normoxic conditions. The ATM inhibitor was kindly provided by Merck KGaA, and dissolved in DMSO (PAN-Biotech, Aidenbach, Germany) and diluted in RPMI 1640 medium. The inhibitor was added to H460 and H1437 cells at 100, 200, or 500 nM just before incubation under hypoxia or normoxia started. Controls also contained < 0.1% DMSO. Inhibitors were left in the media for 24 h and then replaced with fresh RPMI 1640 medium and the plates were returned to the incubator for colony formation. After 4 days (A549), 5 days (H460, Renca, and B16), or 7 days (H1437), plates were imaged by an online microscopy system at 4x magnification (IncuCyte, Essen Bioscience, Sartorius, Göttingen, Germany). The images were analyzed by the IncuCyte Zoom Software (ver. 2016a) (Essen Bioscience, Sartorius, Göttingen,) and colony counts were confirmed by manual curation.

Dose Planning and Simulations
The irradiation plan was carried out as a step and shoot intensity-modulated radiotherapy (IMRT) plan, describing the different dose levels as separate target regions. Planning was done with Raystation treatment planning system (RaySearch Laboratories, Stockholm, Sweden) based on a CT scan of the hypoxia chamber containing 96-well plates filled with water. Irradiation was performed on a Siemens Artiste (6 MV) (Siemens, München, Germany).

Modeling Approach
Large parts of the general model and its derivation were presented and discussed in detail in [11]. Computationally, the code was fully rewritten in Python and elements of GPU computation were introduced. In short, for low LET radiation a homogeneous deposition of dose throughout the cell nucleus with a cell line independent DSB induction rate α DSB = 5 · 10 −3 DSB/(Mbp · Gy), constant over the clinical dose range, wass assumed [36][37][38][39]. Thus, the expected number of DSB in the nucleus ( N tDSB ) can be expressed as: where DNAc is the DNA content of a cell in Mbp and D the applied dose in Gy. The total number of giant loops (N gl ) with a DNA content of DNA gl inside the nucleus is then given by: In this work we assumed DNAc and DNA gl to be 6 Gbp and 2 Mbp, respectively. We implemented a Monte Carlo routine in which, at each iteration, the number of total DSB in the nucleus (N tDSB ) was sampled following a Poisson distribution with the expectation value given by Equation (1). After randomly distributing the sampled amount of DSBs over the giant loops, the number of giant loops without any DSB (N 0 ), with an isolated DSB (N iDSB ), or a complex DSB (N cDSB ) were scored. With the lethality parameters K iDSB and K cDSB , which represent the probabilities of an isolated lesion and a complex lesion leading to cell death, respectively, the probability of the cell to survive (S) is given by [13]: We obtained the expected fraction of a cell population surviving an irradiation by meaning S values obtained from the Monte Carlo algorithm. The lethality parameters can be determined by fitting the result of this routine to survival data. Experimental evidence suggests that one can assume a homogeneous distribution of oxygen inside the nucleus, no change in DNA content or nucleus volume under reduction of oxygen supply [40], and that the oxygen concentration in the cell does effect the initially induced total number of DSB but not their repair rate [6]. Thus, in our model, a change in oxygenation solely leads to an introduction of a modified DSB induction rate α O 2 DSB , which is given by: where α DSB is the rate under normoxia. The modification of α DSB leads through Equation (1) to a change of N tDSB to N O 2 tDSB , which subsequently leads to alterations of N iDSB and N cDSB to N O 2 iDSB and N O 2 cDSB , respectively. Keeping the lethality parameters invariant, as implied above, the survival probability of a cell under hypoxic conditions can then finally be written as: The HRF O 2 DSB value can be determined by fitting the model to hypoxic data, while keeping K iDSB and K cDSB at the values derived from normoxic data. If either one, hypoxic or normoxic data, is not available, the HRF O 2 DSB can be estimated from the parameterization: introduced in our previous work [11], proposed by Carlson et al. [6] and inspired by the initial works of Alper and Howard-Flanders [28]. Fitting this parameterization to data available in literature, we found the values m = 2.94 and K = 0.129% [11]. In order to model the increased cell killing of repair-deficient mutant cell lines or cells exposed to different concentrations of a repair inhibiting drug, a radiosensitization factor (RSF) was introduced into the model. The RSF modifies the lethality parameter of isolated damages K iDSB only, so that the survival probability of a repair impaired cell reads: The RSF is introduced as a free parameter, which is determined by fitting the result of the modified survival probability given in Equation (7) to repair impaired data, while the lethality parameters K iDSB and K cDSB remain set to the values derived from wild type/non-treated cells. As argued above, no interaction between the oxygenation status and the repair capacity is assumed. Therefore, also in the case of a modification of the isolated lesions lethality, it is assumed that the RSF value is set to be invariant under change of the oxygen concentration.