A Multi-Compartment Model of Glioma Response to Fractionated Radiation Therapy Parameterized via Time-Resolved Microscopy Data

Purpose Conventional radiobiology models, including the linear-quadratic model, do not explicitly account for the temporal effects of radiation, thereby making it difficult to make time-resolved predictions of tumor response to fractionated radiation. To overcome this limitation, we propose and validate an experimental-computational approach that predicts the changes in cell number over time in response to fractionated radiation. Methods We irradiated 9L and C6 glioma cells with six different fractionation schemes yielding a total dose of either 16 Gy or 20 Gy, and then observed their response via time-resolved microscopy. Phase-contrast images and Cytotox Red images (to label dead cells) were collected every 4 to 6 hours up to 330 hours post-radiation. Using 75% of the total data (i.e., 262 9L curves and 211 C6 curves), we calibrated a two-species model describing proliferative and senescent cells. We then applied the calibrated parameters to a validation dataset (the remaining 25% of the data, i.e., 91 9L curves and 74 C6 curves) to predict radiation response. Model predictions were compared to the microscopy measurements using the Pearson correlation coefficient (PCC) and the concordance correlation coefficient (CCC). Results For the 9L cells, we observed PCCs and CCCs between the model predictions and validation data of (mean ± standard error) 0.96 ± 0.007 and 0.88 ± 0.013, respectively, across all fractionation schemes. For the C6 cells, we observed PCCs and CCCs between model predictions and the validation data were 0.89 ± 0.008 and 0.75 ± 0.017, respectively, across all fractionation schemes. Conclusion By proposing a time-resolved mathematical model of fractionated radiation response that can be experimentally verified in vitro, this study is the first to establish a framework for quantitative characterization and prediction of the dynamic radiobiological response of 9L and C6 gliomas to fractionated radiotherapy.


INTRODUCTION
Radiation therapy is a central component of the standard-of-care for treating malignant gliomas (1), especially when the tumor is located near sensitive brain regions with important functions that are unresectable by surgery. Though various dose escalation and fractionation schemes (i.e., hyper-and hypo-fractionation) have been investigated, none have shown definitive improvement on the long-term survival for glioblastoma patients (2). One reason for this limitation is that the efficacy of radiation therapy varies between patients due to heterogeneous radiosensitivity of the cells within each individual's tumor (3). If there was a mathematical model that could accurately characterize, and predict, the response of tumor cells to radiation therapy with patient-specific data, then there would be the opportunity to optimize the radiation plan for each individual (4). The currently accepted model for evaluating radiation response given a specific dose is the linear quadratic (LQ) model which was originally developed empirically more than 40 years ago (5). The LQ model quantifies the survival fractions of cell colonies given a specific radiation dose and, though it provides a simple, and practical relationship between those two measureables, it is not without its limitations. In particular, the LQ model does not explicitly characterize the temporal changes in tumor cell number; that is the LQ model is not a function of time. Thus, while it can provide accurate predictions of endpoint predictions (6), it is not capable of predicting the temporal dynamics of radiation response. Additionally, interpretation of the two main parameters in the LQ model (alpha and beta) is fraught with difficult, thereby clouding their biological meaning (7). This is despite the now vast biological knowledge that exists regarding DNA repair (8) and radiation-induced cell death pathways (9). To address these two limitations, we previously proposed and validated a mechanism-based time-resolved model (10) to a single-dose treatment. We now seek to extend this model to account for multiple-fraction treatment regimens.
Though a large single-dose of radiation can effectively kill tumor cells, it is rarely used in clinical settings as high doses also cause irreversible cytotoxicity to surrounding healthy tissues. Thus, the notion of delivering a target total dose in "fractions" over an extended period of time was adopted. There are four key conceptions that are frequently kept in mind when designing multi-fraction treatment plans: DNA damage repair, repopulation, cell cycle redistribution, reoxygenation [sometimes referred to as the "Four R"s (11)]. The DNA damage repair mechanisms help the nearby healthy tissue recover between treatment intervals (12), with the hope that the repair mechanisms are erroneous in the tumor cells leading to their eventual cell death after repeated fractions (13). Tumor cell repopulation (i.e., the ability of tumor cells to proliferate between treatment intervals) can undermine radiation efficacy, and thus may require extra fraction and/or total dose to achieve tumor control (14). Cell cycle redistribution increases the average tumor cell killing by allowing radiationresistant cells in S phase to redistribute into the more sensitive M phase (15). Reoxygenation also enhances radiation damage as radiation can produce free radicals which damage DNA, and this damage can be made permanent by the presence of molecular oxygen [i.e., the 'oxygen-fixation hypothesis' (16)]. Additionally, hypoxic regions of a tumor regions may become reoxygenated between fractions (17). Previous modeling work has focused on quantifying the effects of these four "R"s on the endpoint survival fraction by (for example) incorporating an "oxygen enhancement ratio" (18) or repopulation (19) into the LQ model. More recently, several studies have constructed radiation response models that account for temporal changes in hypoxia (20), DNA repair (21), and fractionation (22). Hormuth et al. contributed a tissue-scale model that employed the oxygen enhancement ratio to adjust the radiation efficacy during fractionation treatment and tested model predictions against in vivo MRI data (23,24). Brüningk et al. (25) proposed a cell-scale decision tree model that accounted for conversion between cell cycle compartments after radiation. All these models indicate an increasing interest in mathematically describing the temporal dynamics of radiation response that are not captured by the conventional LQ-based models.
We first propose to extend our previous single-dose model (10) to characterize the radiation response of gliomas to fractionated treatment. The fractionation model explicitly incorporates temporal changes due to DNA damage repair, cell repopulation, and cell cycle effect related to senescence. We then perform in vitro microscopy experiments with 9L and C6 cell lines to obtain the radiation response curves collected at high temporal resolution under different treatment schedules and total radiation doses. Our model is then trained on 75% of the total data to calibrate the parameters. Finally, the remaining 25% of data serve as a validation group to assess the model's predictive accuracy. Our mechanismbased, time-resolved, mathematical model achieves high predictive accuracy across a range of fractionation schedules verified by six different fractionation schemes and both cell lines.

Cell Culture
The 9L (American Type Culture Collection, ATCC) and C6 (Sigma Aldrich) cells are cultured according to the manufacturer's guidelines as previously described (10). The 9L cell line is cultured with Eagle's minimum essential medium (ATCC, VA), and the C6 cell line is cultured with Ham's F12 (Corning, NY). Both cell lines media are supplemented with 10% FBS and 2 mM L-Glutamine. 0.2% Plasmocin Prophylactic (In vivogen, CA) is supplied in the media to prevent mycoplasma contamination. The mean passage number of the cells used in the experiments was 50 ± 20. Both 9L and C6 are rat cell lines, but are commonly used for general glioma studies (26, 27). Figure 1 illustrates the radiation treatment schedule employed in these studies. 9L and C6 cells were seeded on 96-well plates (Corning, NY) at densities ranging from 3,200 to 32,000 cells/cm 2 (1,000 to 10,000 cells total per well). To avoid cells reaching the carrying capacity at later timepoints (which can result in cell death due to lack of nutrients and physical space), we do not seed at a confluence higher than 10,000 total cells. The cells are then incubated overnight (~12 hours) to allow for attachment and recovery. Before irradiation, the media is changed and augmented with 250 nM Cytotox Red dye (Cat. No. 4632, Essen BioScience, MI), a non-perturbing fluorescent dye to label dead cells. For both the 9L and C6 cell lines, we separate the wells into either a 16 Gy or a 20 Gy total dose group. In the 16 Gy group, we irradiate the cells with either four fractions of 4 Gy, three fractions of 5.3 Gy, or two fractions of 8 Gy with 24-hour intervals between every fraction. In the 20 Gy group, we irradiate cells with four fractions of 5 Gy, three fractions of 6.7 Gy, and two fractions of 10 Gy with 24-hour intervals between every fraction. All radiation is delivered by a CellRad irradiator (Faxitron X-Ray Corp, Wheeling, IL, MA) at a dose rate of 1.5 Gy/min (130 KeV, 5 mA, 0.5 mm aluminum filter). After treatment, phase-contrast images and fluorescent Cytotox Red images (for labeling dead cells) are acquired immediately after the first fraction via the Incucyte S3 live imaging system (Essen BioScience, Ann Arbor, MI) with a 4× objective, whole-well imaging mode every four to six hours up to approximately 330 hours post-irradiation. Our media culture, along with the Cytotox Red dye, was refreshed every five days throughout the experiment. To prevent cell loss when refreshing the media in the 96-well plates, we pipetted only the top 80 ml of the total 100 ml per well to minimize the disturbance to attached cells. Live and dead cells were segmented using a semi-automated pipeline consisting of a histogram Otsubased method followed by a morphology-based cell debris removal [described in detail in (10)].

DSB Repair Kinetics
To quantify radiation-induced DNA double strand breaks (DSB) and repair kinetics, we previously measured the expression level of the gH2AX protein [a commonly used DSB biomarker (28)] via flow cytometry after irradiating cells with a single dose of 2, 4, 8, or 16 Gy (see the Supplementary Materials of (10) for details). We then used linear interpolation to obtain the DSB repair kinetics for all other doses. The same data is used in this study as described below.

DNA Repair
DNA repair is represented by an exponential decay equation: where f DSB (t,D) is the fraction of DSBs remaining unrepaired (normalized between 0 and 1) at time, t, k repair (D) (in units of hr -1 ) is the DSB repair rate at dose per fraction D. We use the same k repair measured from our single-dose treatment study as an estimate of the DSB repair rate during fractionation schedules. This assignment is justified as Mariotti et al. (29) have measured gH2AX under both single and multi-fractionated treatment schemes and showed similar repair kinetics in response to a second dose when cells are given proper time for repair and recovery. Our previous flow cytometry experiments (10) indicate that most (> 80%) DSBs are repaired within 24 hours given the dose range we employ in the experiments. Therefore, this estimation is reasonable.

Mathematical Modeling of Cell Growth and Response to Radiation Therapy
The fractionated treatment model is an extension of our previous single-dose radiation model described in (10) that models cell response as a function of early cell death (corresponding to apoptosis) and late cell death (corresponding to mitotic catastrophe). We present the salient details of our single-dose radiation model, but note that the complete development and underlying assumptions are detailed in (10).

Single Species Model of Cell Growth in the Absence of Radiation Therapy
For glioma cell proliferation in the absence of radiation therapy, we augment exponential growth by incorporation of the logistic growth and Allee effect: |{z} logistic growth (2) where N(t) is the tumor cell confluence, k p is the proliferation rate (see Table 1 for a listing of all model parameters, their definition, and units), A quantifies the strength of the Allee effect, and q is the carrying capacity (i.e., the maximum number of cells that can fit within a given volume due to space and nutrient limitations). The Allee effect describes the cooperation effects in cell proliferation rate and is proved significant in glioblastoma progression by Neufeld et al. (30). Our previous analysis (10) used model selection to establish that both logistic growth and the Allee effect are necessary to accurately describe our glioma data.

Single-Species Model of Cell Growth in the Response to Radiation Therapy
After radiation, a small number of cells undergo early apoptosis, which is an outcome of activation of DNA protein kinase and p53 due to excessive DSBs (31); these events occur on the timescale of hours to days (32). Meanwhile, DNA misrepair does not directly kill cells and can accumulate within cells' genome, eventually triggering chromosome aberration and mitotic catastrophe (33,34); these events occur on the timescale of days to weeks (32). Based on the above mechanisms, we add early and late death terms to Eq. (2): where k ed and k ld (both in units of hr -1 ) represent early death and late death rates, respectively, and are a function of time t, dose per fraction D, and the initial confluency N 0 . We use the following equation to describe early cell death as a function of the fraction of unrepaired DSBs within the cell's genome, f DSB (t): where k acute (D,N 0 ) is the acute death rate. t i (units: hours) is the time that has passed since the i th fraction. As it would be extremely complex to model the death due to each fraction separately, k acute (D,N 0 ) is viewed as an average death rate across all fractions and written as: where a acute , N is a scale factor representing the contribution of acute death due to the initial seeding density, N 0 . Biologically, N 0 influences (due to cell-cell contact) the proportion of actively proliferating cells at the time of radiation, which determines radiation sensitivity[(as late G2 and M phase is the most sensitive cell cycle (35)]. k acute , D is a death rate indicating the contribution of radiation dose to acute death, as larger doses can translate to a higher number of DSBs and, therefore, early apoptosis. Figure 2 illustrates the changes in k ed over time. We note that the "+1" in Eq. (6) is for mathematical convenience. When a acute , N (i.e., the contribution of the initial seeding density to acute death) is close to 0, Eq. (6) simplifies to k acute (D,N 0 )= k acute , D ·D. Thus, without the "+1", when a acute , N approaches 0, k acute (D,N 0 ) will also decrease to 0 and the effect of dose will only be determined by seeding density effects. a accum,N 1 Scale factor quantifying the contribution of initial confluence to late death k accum,D hr -1 Death rate quantifying the contribution of radiation doses to late death r hr -1 Radiation efficacy k ps hr -1 Conversion rate from proliferation to senescent components.
The unit "1" means the parameter is unitless.
The "late death" component models mitotic catastrophe of misrepaired cells that occurs following several cell divisions. This population has previously also been referred to as the "abortive fraction" (36). Some cells can survive hours to weeks after radiation before mitotic catastrophe occurs. Here, we model late death as: where k accum (D,N 0 ) is that rate at which the misrepair error first accumulates within the cells' genome, and r (in units of hr -1 ) controls the decay rate of the radiation efficacy. The decay rate of radiation efficacy is viewed as an intrinsic property of each cell line and is a constant across different treatment conditions. Just as for the early death term, we model the parameter k accum (D,N 0 ) as an average across fractions, instead of modeling each fraction separately. The accumulation death rate k accum (D,N 0 ) is thus written as: where a accum , N is a scale factor representing the contribution of late accumulation cell death due to the initial seeding density, N 0 . Eq. (8) accounts for the fact that cell density directly impacts the proliferation rate via (for example) the proportion cells that are in M phase, which eventually determines how many cells can undergo mitotic catastrophe as it occurs during M phase. k accum,D is a death rate describing the contribution of radiation dose to accumulation death, as high doses induce a high probability of misrepair [caused by misjoining of clustered DSBs (37)]. Note that Eqs. (6) -(8) have a modified form from our previous single dose study (10); we return to this point in Discussion section.

Two-Species Model of Cell Growth in Response to Radiation Therapy
We construct a two-species model of response to radiation therapy by considering both proliferative and senescent tumor cells. Several mechanisms in radiobiology contribute to the appearance of a senescent population after radiation (38) including (for example) cell cycle checkpoint pathways (39). These senescent cells can remain metabolically active, but undergo irreversible cell cycle arrest and thus can no longer replicate. Without this component, the model assumes cells either return to proliferating (overestimate cell survival) or undergo early or late death (overestimate radiation cell killing), thereby causing a systematic error in the predicted confluence. Therefore, we extended Eq. (3) to include a proliferating compartment, N p (t), and a senescent compartment, N s (t): Eq. (9) -(10) characterize the conversion from the proliferation to senescent compartment due to radiation. As cell-cell contact promotes senescence, we assume that the convertion rate follows a simple linear relation proportional to the initial cell density with the proportionality constant k ps (units of hr -1 ). Note that the carrying capacity, q, is now shared by both N p and N s . Additionally, the conversion rate k ps is a function of the total dose (i.e., 16 Gy or 20 Gy in our experiments). This simplification (i.e., making k ps a constant based on the total dose) is because we do not have a direct measurement of the senescent population as this population is changing with each dosing scheme and over time.
That is, k ps should really be a function of time, dose per fraction, and fraction number. To practically realize such an explicit expression for k ps would require additional experimental measurements.

Numerical Implementation of the Mathematical Models
ODE models were implemented via the finite difference method with a fully explicit forward Euler formulation with a time step of 0.01 hrs with initial condition, N p (0), equal to the cell confluence measurements at time 0 and N s (0) = 0. Early acute death rate term is multiplied by a smooth heaviside step function, via the hyperbolic tangent function; tanh [f DSB (t)], to prevent curve discontinuity.

Model Selection
Eqs. (1) -(10) are based on general radiobiology mechanisms. For specific cell lines with different signaling pathways and radiation sensitivity, it is unclear if this model is appropriate to describe the data. Therefore, it is crucial to perform model selection prior to applying the model to a specific cell line. Starting from the above "full model", we systematically remove one or two parameters, yielding seven competing "daughter" models. Using the early death term as an example, the initial seeding density N 0 and doses D might not have an impact on early apoptosis for the 9L and C6 cell lines. By removing either a acute , N or k accum,D or both, we obtain three "daughter" models (i.e., models 2-4 in section 2 of the Supplementary Materials).
Similarly, removing the late death term (three models) or the senescent term (one model) yields an additional seven "daughter" models. (See the section 2 of the Supplementary Materials for the formulation of all eight models.). For each model, parameters are fit with the Levenberg-Marquardt algorithm (via "lsqnonlin" in MATLAB). To determine which mechanisms are required for optimally characterizing 9L and C6 data and to obtain the most parsimonious model, we perform model selection on these eight models via the Akaike information criterion (AIC) (40): where n is the sample number (i.e., the number of cell confluency curves in our scenario), RSS is the residual squared sum between the model fit and data, and p is the number of free parameters. The AIC finds the most parsimonious model by balancing the relative goodness-of-fit with the number of free parameters. Specifically, we build our training set by randomly picking seventy-five percent of the data under each treatment condition, leading to 262 replicates for the 9L cells and 211 replicates for the C6 cells. (That is, we pick 75% from four fractions of 4 Gy, 75% from two fractions of 10 Gy, etc., thereby ensuring that each treatment condition is equally represented in the training set.) The remaining 25% of the data is used for validation. (Note that the word 'replicate' is in reference to our experiments; that is, we repeated a group of independent wells with the same treatment schedules at time point 0. The response from each well over the 340 hours won't be identical as they are affected by (for example) subtle variations in the cells' phenotype or genotype, variations in initial seeding density, etc.) The AIC-based model selection is then performed on all eight models including the full model using the training set. By globally fitting the training set (i.e., 262 9L replicates and 211 C6 replicates, to these eight models respectively) we obtain the residual squared sum over the entire training set for each model. The model that returns the lowest AIC score is selected as the most parsimonious. We also compute the Akaike weights via: where AIC i is the AIC score of the i th model, and AIC min is the minimum AIC observed among all eight models. These weights are used to compare the models to each other. Table 1 lists parameter definitions and the methods we use to obtain these parameters. We previously (10) fit untreated cell data to Eq. (1) to obtain the proliferation rate, k p , carrying capacity, q, and the Allee constant, A, for both the 9L and C6 cells. These parameters are assumed constant throughout the present study. The standard deviation of the estimated model parameters is computed via "nlparci" in MATLAB, and employed to generate a corresponding distribution via "makedist". As the biological definitions of these parameters specify their values must be positive, the parameters are given a lower bound of zero during calibration. Note that during model calibration, we fit parameters globally in the sense that all curves from the training set, consisting of both the 16 Gy and 20 Gy total dose groups, are fit together since all parameters are independent of dose and initial cell density except for the conversion rate k ps . As k ps is a function of the total dose, we treat k ps at 16 Gy and at 20 Gy as two individual parameters in our calibration process. Distributions of the calibrated parameters are compared between the two cell lines to verify if they are significantly different by the z-test (via "ztest" in MATLAB at the 5% significance level).

Model Validation and Error Analysis
Validation is performed on the AIC selected model using the remaining twenty-five percent of data (i.e., 91 9L curves and 74 C6 curves), which are "unseen" during both the model selection and parameter calibration steps. The forward model has two inputs: initial seeding density (N 0 ) and dose schedule. We run the model forward using the calibrated parameters as shown in Figure 5. The prediction mean and intervals are then computed via the MATLAB function "nlpredci" by inputting the calibrated parameters, Jacobian matrix, and residuals determined by the "lsqnonlin" during calibration. Radiation responses are predicted using the same set of parameters regardless of treatment dose schedules or initial density. For example, predicting the effects of the four fractions of 4 Gy on the low confluence group employs the same set of parameters as predicting the effects of the two fractions of 10 Gy on the high confluence group (except for the k ps , which is a constant based on total dose). To quantify the predictive accuracy, we compare the measured and model prediction means using the Pearson correlation coefficient (PCC) and concordance correlation coefficient (CCC).

Image Segmentation and Cell Response Curves
The same image segmentation pipeline from our previous study (10) is used here. When comparing the automatically segmented images to our manually segmented baseline, this segmentation pipeline achieves an average Sørensen-Dice coefficient of 0.79 (i.e., 21% error). See Figure 3 for an example of a cell response curve that received four fractions of 4 Gy and its corresponding image segmentation. Figure 4 summarizes the AIC weights across all eight models. As it has the highest weights for both the 9L and C6 cell lines, Model 3 was selected as the most parsimonious model from the training set and thus will be used for parameter calibration and  prediction. While the complete formulation is presented in the section 2 of the Supplementary Materials, Model 3 combines early apoptosis (k acute,N ), mitotic catastrophe (a accum,N , k accum,D , r), and senescence [k ps (16 Gy), k ps (20 Gy)] as given by the following system of equations:

Model Selection
Note that the parameter k acute , D in Eq. (6) is removed and the corresponding expression for early death is rewritten as in Eq. (15). This was done since model selection indicated that k acute , D is not required to characterize our data. Consequently, we removed k acute , D and modified the notation of a acute , N in Eq. (6) to k acute , N (in units of hr -1 ) in Eq. (15) to represent the death rate. From the AIC score, model 3 [i.e., Eqs. (14) - (17)] is 1.18 times more likely to be the best model than model 4 'no early death' model, 2.61 times more likely than model 7 'no late death model', and 1.82 times more likely than model 8 'no senescence' single species model. These results indicate the importance of accumulation effects (i.e., accumulating DNA misrepair, mitotic catastrophe and gradual conversion to senescence) over acute effects (i.e., early apoptosis) to quantify the time courses of 9L and C6 radiation response. Note that the goal of the model selection process is to identify the most parsimonious model to describe our time-resolved data, rather than the model that includes the most biology. As for the 9L and C6 cell lines, most cells die due to late mitotic catastrophe while early apoptosis only kills a minority of cells. Thus, removing the radiation dose effect from the early death term as in model 3 does not harm the model's ability to characterize the data.

Parameter Calibration
The AIC selected model has six parameters k acute,N , a accum,N , k accum,D , r, k ps (16 Gy), and k ps (20 Gy). The top panel in Figure 5 shows the calibrated parameters for the 9L cells, while the bottom panel shows the calibrated parameters for the C6 cells. Note that the model allows to quantify the degree to which the C6 cell line is more radiation sensitive than the 9L cell line, as has been previously reported (41). In particular, the early death rates k acute,N (p value = 9.5e-23), late death rates k accum,D (p value = 8.1e-21) and conversion rate to senescent component k ps (p value = 1.6e-20 for 16 Gy and 0.0 for 20 Gy) of the C6 are all significantly larger than 9L via the z-test. Both the 9L and C6 exhibit a similar radiation efficacy decay r (p value = 0.44, i.e., no significant difference via z-test), suggesting a similar duration of radiation cell killing persisting on both cell lines. The proliferation rate, carrying capacity, and Allee effect parameter values, as well as the number of training curves for calibration, are provided in the section 1 of the Supplementary Material.

Model Validation and Error Analysis
We use the validation group (25% of our total data, 91 9L curves and 74 C6 curves) to evaluate the predictive accuracy of our model. Figure 6 presents examples of the model validation results with each column representing one initial confluence of a specific cell line (e.g., the first column shows a low initial seeding confluence of the 9L cell line). Each row shows a different fractionation schedule (e.g., the first row shows cells receiving four fractions of 4 Gy radiation). The error bar on the measurement (labeled by blue) is based on the image segmentation error (21%) from our previous study (10), as the same segmentation pipeline is used. The prediction error (labeled by red) is computed from MATLAB function 'nlpredci', which computes the prediction interval via the Delta method based on the Jacobian matrix. For the 9L group receiving a total dose of 16 Gy, the average PCC and CCC between the predicted and the measured data are 0.92 ± 0.009 (average ± standard error) and 0.78 ± 0.014, respectively. For the 9L group receiving a total of 20 Gy, the PCC and CCC are 0.98 ± 0.001 and 0.96 ± 0.002, respectively. For the C6 group receiving a total of 16 Gy, the PCC and CCC are 0.89 ± 0.011 and 0.77 ± 0.020, respectively. Finally, for the C6 group receiving a total of 20 Gy, the PCC and CCC are 0.90 ± 0.004 and 0.73 ± 0.014, respectively. Details of PCC and CCC values for each treatment condition is provided in Table 2.
Overall, the accuracy of prediction is superior for the 9L cell line than the C6 cell line. This is at least partially due to the heterogeneous radiation response observed across the C6 replicates; we return to this important point in Discussion section. FIGURE 5 | Parameter calibration results. All parameters are fit globally using the training set and are independent of initial seeding density or dose schedules (e.g., the four fractions of 4 Gy curves share the same set of parameters as the two fractions of 10 Gy curves for both the 9L and C6; the one exception is for k ps , where we treat k ps (16 Gy) and k ps (20 Gy) as two individual parameters. The biological interpretation of the parameters requires that these parameters must take on a value great than or equal to 0; thus, the lower bound of the parameters is set to zero during calibration. FIGURE 6 | Model validation. The prediction interval (labeled by red) and measured microscopy data (labeled blue dots) are plot for representative examples of each initial and treatment condition. Previously we determined our average segmentation error was 21% using the Sørensen-Dice coefficient (10). This segmentation error is indicated by the blue intervals. Each row in this figure shows different treatment schedules labeled on the left; for example, the first row shows cells treated with four fractions of 4 Gy. Each column stands for the initial confluence for a specific cell line; for example, the first column represents low initial seeding density for 9L cell line. Predictions are made globally for each of the cell lines; that is, predictions for the 9L cells with an initial low confluence receiving four fractions of 4 Gy are made using the same set of parameters as the 9L cells with an initial high confluence receiving two fractions of 10 Gy. Given the initial confluence and treatment schedule as the two inputs, our model makes accurate predictions across a wide range of initial conditions. However, the model is not perfect as the prediction typically undershoots the first peak while overshooting the tail for C6 cell line; an important point we return to in the Discussion section.

DISCUSSION
We have proposed an experimental-computational system that includes a mathematical model that characterizes the dynamic cell response of the 9L and C6 glioma receiving fractionated radiation. Model parameters are calibrated using a training set of various initial seeding densities and dose schedules over the course of two weeks. Next, the calibrated parameters are applied to a validation set to predict response to radiation therapy. Each term in our model has an explicit or implicit (i.e., a accum,N , k ps (16 Gy), k ps (20 Gy) are indirectly related to cell cycles) biological definition. For example, early apoptosis and late mitotic catastrophe are captured by the parameters k acute,N and k accum,D , r, respectively. This approach has several advantages over current models as it explicitly includes the temporal dynamics of several biological mechanisms related to the response of cells to fractionated radiation.
Mathematically, the LQ model describes the fraction of survival cells at experimental endpoints given a specific radiation dose, and therefore is not designed for calibrating with time-resolved data. Though there are mathematical models that account for temporal dynamics by embedding the LQ term into the death rates (42), the true death rates are unlikely to obey the LQ relation throughout the whole time course as it ignores time dependent phenomena such as early apoptosis and late mitotic catastrophe. Timedependent radiobiologic mechanisms are often ignored in these modeling studies, a limitation possibly due to the historical difficulty of accessing radiation response data with high temporal resolution. However, recent advances in imaging techniques can now provide the requisite, high temporal resolution data (43)(44)(45) appropriate for model calibration. When we mathematically characterize these time series data, a dynamic model is required, since the measurement interval (e.g., hours) is much shorter than the length of the experiment when the surviving fraction would be determined via the LQ model. To the best of our knowledge, this study is the first mechanism-based fractionated radiation model verified by the cell experimental data at high temporal-resolution. We have constructed this approach by building on previous time-resolved dynamic radiation models that either calibrated to a single dose of radiation (46), or to data collected every several days (23) or weeks (47). We hope that the current effort can provide motivation for focusing on factors that describe the dynamics of radiobiology as a function of time, thereby potentially providing a new way to guide and optimize radiation dose scheduling.
As each term in the selected model system [i.e., Eqs. (14) - (17)] is based on an underlying mechanism, the calibrated parameters summarized in Figure 5 are easily interpreted in light of the underlying biology. Though the specific values of parameters may differ for each cell line (a not unexpected observation), the trend observed (i.e., radiosensitive cells exhibit significantly larger death rates) may carry over to other cell lines since the proposed model and the parameters are based on biological mechanisms common to a wide range of cell lines (48). Thus, the model summarized by Eqs. (14) - (17) is likely applicable to other cell lines for which radiation treatment is of interest, thereby having a significant impact beyond gliomas.
Despite the advantages, there are several opportunities for improvement in both the experimental and modeling components of the study. For example, note that in Figure 6 some curves suffer a discontinuous jump (e.g., at approximately 110 hours in the "4 fractions of 5 Gy" high confluence group). This is due to the cell loss when we refresh the media. While we have implemented a method to change the media as delicately as possible (as described in section 2.1.2), we are still aspirating an unknown number of live cells at each media change. Potentially more significant, though, is how we handle the senescent component. Since we did not have a method in place that allows us to directly measure this component over time, we are forced to infer its existence and behavior via parameter calibration. Clearly, incorporating longitudinal measurements of the fractions of senescent cells is required to generalize the model to other dosing schemes.
Areas for improvement on the modeling side include developing a more rigorous linking between radiation dose and the rate constants. For example, note that Eqs. (6) and (8) have a different form from our previous single-dose study (10). The main reason for the changes is due to the dose range in this study (4 -10 Gy per fraction) compared to the previous study (a single fraction of 2 -16 Gy). In our previous single-fraction study we observed a saturation effect in the death rates; i.e., the death rates do not significantly increase above a certain dose threshold. This is because cell death does not happen instantaneously; i.e., the cell death pathways require time to be executed. Thus, we use Michaelis-Menten equations to characterize the observed saturation effect. However, in the present study, we use a narrower range of radiation dose and we did not observe the saturation in cell killing within this dose range. Thus, we use a simple linear relation for the acute death, k acute (D,N 0 ), accumulation death, k accum (D,N 0 ), and a constant radiation efficacy r, versus radiation doses. Consequently, the death rates k acute (D,N 0 ), k accum (D,N 0 ), and r most likely require future modifications when applied to dose schedules involving larger doses. We also note that, compared to the previous single-dose study [i.e., Eq. (6) in (10)], we removed the parameter T in Eq. (7). This parameter assumes there is a time delay, T, before the accumulation of misrepairs can start. Our previous study established its value between the first 20-40 hours after radiation for the 9L cell line, and 0 hours for the C6 cell line [consistent with the established radiosensitivity of the C6 cell line (41)]. In current fractionated dosing scheme, as we extend the experiment time to approximately 340 hours, the value of T is much less than the length of the experiment. Thus, it is reasonable to remove the parameter T for simplicity. Another limitation is that Eqs. (14) - (17) do not (of course) account for all of the most prominent radiobiological processes. For example, we assume a fixed conversion rate between the proliferative and senescent components-a parameter that is likely to change with time, dose per fraction, and fraction numbers, and therefore would require a function of its own to describe its temporal dynamics. Indeed, this may explain the undershoot peak or overshoot tail in Figure 6. A second area for improvement in Eqs. (14) - (17) is the explicit incorporation of phenotypic or genomic heterogeneity across the population, which almost certainly affects the early (k acute,N ) and late (k accum,D ) death effects between different replicates. This explains why the correlation coefficients in Table 2 indicate the model performs worse under certain dose schedules for the C6 cell line. See section 3 of the Supplementary Materials for an example of the heterogeneous response we observe in two replicates of C6 cells treated with the same three fractions of 6.7 Gy (note this is the group with the lowest PCC and CCC value). Accounting for such heterogeneity to improve the predictive accuracy of the model will be the focus of future study.
As presented, the current model has limited clinical application because it is formulated for describing the 2D dynamics of cells in a dish which is, of course, very different than the in vivo situation. However, the concept of using time-resolved data to calibrate a biologically-based, mathematical model to make patient specific predictions is highly translatable and, indeed, something we have investigated at length in both the in vivo pre-clinical (49)(50)(51)(52)(53) and clinical (24,(54)(55)(56)(57)(58) settings.

CONCLUSION
We have extended our previous single-dose model to account fractionated treatment, and successfully validated the resulting model with in vitro experimental microscopy data. This study demonstrates a promising experimental-mathematical approach based on radiobiology mechanisms that can accurately predict the temporal dynamics of the response of glioma cells to radiation. Future efforts include linking the model to our tissue scale formalism for predicting response in patients (23), and employing the methods of optimal control theory to optimize treatment outcomes in pre-clinical murine studies (59).

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
JL collected data, performed modeling analysis and wrote the manuscript. JY helped with the experiments. DH and TY conceptualized and supervised the study, reviewed and edited the manuscript. All authors have read and approved the final manuscript.