Abstract
In this work, a method is established to calibrate a model that describes the basic dynamics of DNA damage and repair. The model can be used to extend planning for radiotherapy and hyperthermia in order to include the biological effects. In contrast to “syntactic” models (e.g., describing molecular kinetics), the model used here describes radiobiological semantics, resulting in a more powerful model but also in a far more challenging calibration. Model calibration is attempted from clonogenic assay data (doses of 0–6 Gy) and from time-resolved comet assay data obtained within 6 h after irradiation with 6 Gy. It is demonstrated that either of those two sources of information alone is insufficient for successful model calibration, and that both sources of information combined in a holistic approach are necessary to find viable model parameters. Approximate Bayesian computation (ABC) with simulated annealing is used for parameter search, revealing two aspects that are beneficial to resolving the calibration problem: (1) assessing posterior parameter distributions instead of point-estimates and (2) combining calibration runs from different assays by joining posterior distributions instead of running a single calibration run with a combined, computationally very expensive objective function.
1. Introduction
DNA damage and repair is a critical aspect of radiotherapy, where tumor cells are killed by irradiation. The radiation induces DNA damage which eventually leads to cell death if the damage cannot be repaired successfully. Mild hyperthermia is a treatment to boost radiotherapy by heating up the cancer cells to temperatures between 41°C and 43°C. While the exact working principles of hyperthermia and its interaction with radiotherapy is still subject to research [1], it has been shown that hyperthermia acts as a radiosensitizer by affecting the DNA repair that takes place after an irradiation event [2–5]. In consequence, knowledge about the dynamics of DNA damage and repair is essential in order to optimize hyperthermia treatment plans.
For radiotherapy, in silico modeling is employed to assist in treatment planning decisions. Such planning is based on Monte Carlo simulations or kernel methods and deliver dose-volume histograms [6]. Beyond these geometric dose calculations, approaches to shape the prescribed radiation dose according to the biological properties of the tumor have been proposed but are currently not established [7]. The prescribed dose of radiation is generally divided into fractions that are delivered in subsequent sessions; however, this fractionation scheme is usually not optimized on a patient level, and the dose prescription is chosen based on clinical trials and experience. While planning software may include calculators for biological effective dose (BED) and equivalent dose (EQD2), they are not modeling biological effects (such as DNA damage and repair), but rather, they are tools for comparing fractionation schemes. Similarly, for hyperthermia, planning systems for hyperthermia output temperature or specific absorption rate (SAR) maps exist [8] and calculators for equivalent doses have been proposed [9, 10]. Yet, a more profound understanding and modeling of the aforementioned radiobiological effects—DNA damage and repair in this context—would yield a better treatment method. For example, hyperthermia is believed to deactivate DNA repair proteins for a certain amount of time [2]. If radiation-induced damage is introduced during this time window, odds of eliminating the cells increase [11]. Thus, if calibrated correctly, a model involving DNA damage and repair would be able to quantify the duration of this window by simulating the de- and reactivation of said proteins.
In this work, a method is established to calibrate a model that describes the basic dynamics of DNA damage and repair. This model can then be used to extend planning for radiotherapy and hyperthermia to include the biological effects discussed above, i.e., DNA damage and repair: The biological system is modeled in silico, and a parameter search for model calibration is performed with the goal to be able to quantify biological effects for the system of interest. While previous efforts demonstrated feasibility [12], a thorough analysis of the calibration process is provided. This analysis reveals that some parameters remain unidentified. One strength of the method is that it is able to combine calibration results originating from different input data sources (i.e., assays). With this approach, the yet unidentified parameters could be further refined.
Model calibration requires data which can be obtained from number sources which are shown in Figure 1: (1) immunocytochemical assays such as H2AX, which quantify DNA repair [13]; (2) comet assay, which quantifies the amount of DNA damage [14]; this assay is further discussed in Section 2.3; (3) clonogenic assay, which quantifies clonogenicity [15] and is discussed in Section 2.2. (4) In a clinical setting, DNA damage and repair in tumor cells also affects response evaluation criteria, tumor volume, patient survival, tumor progression and growth rate, etc. Thus, these data (yet quite heterogeneous [7] and thus potentially a poor choice) could, in theory, also be used for model calibration.
These four different options correspond to the four levels illustrated in Figure 1 on the left. On the right, suitable models for these types of readout are depicted. Often, these models merely attempt to replicate some observed readout. For example, the cell survival curves discussed above usually exhibit a parabolic nature in the logarithmic domain [16, 17]. Thus, a quadratic model for is often used for the dose-response, without further rationale but just as a method for fitting the existing data. In the past, this approach has been expanded to a linear-quadratic-linear relationship [18], again in a mere attempt to mimic experimentally observed data. Similarly, biostatistical models for comet assay analysis are able to describe the assay readout but do not model actual DNA damage and repair, let alone in a dynamic fashion [19, 20].
Another class of models go one step further and actually describe underlying molecular principles instead of the mere assay readout. For example, the H2AX phosphorylation discussed above can be modeled using a set of differential equations [21], and the H2AX readout is derived from the model. This approach is mechanistic in the sense that it is directly modeling the kinetics of the H2AX pathway and can be seen as syntactic description of molecular mechanisms. Other models including the lethal-potentially-lethal (LPL) model by Curtis [22] and the model by Vassiliev [23] and the -LQ model [24] all follow radiobiologically motivated approaches but do not consider hyperthermia. The AlphaR model [25] takes the effect of hyperthermia into account, albeit for temperatures above 43.5°C which are not the focus of this work. Going one step further, the multi-hit-repair (MHR) model describes radiobiological semantics [26] instead of mechanics. It was used to derive cell survival curves [12] as well as comet assay readouts [27]. In addition to the semantic approach, the MHR model was chosen for this work because it is bioinspired and in the past, its ability was shown to explain many radiobiological phenomena.
2. Materials and Methods
In the following sections, the experimental setup (Section 2.1), the different biological assays (Section 2.2 and 2.3), the model used in this work (Section 2.4), the methods to map the model state to the readout from experimental assays (Section 2.5), and the calibration method (Section 2.6) are introduced, concluding with a brief section about the software and its availability (Section 2.7).
2.1. Experimental Setup
Hyperthermia and irradiation was performed on cells from the Abrams cell line; they were a kind gift of Prof. Robert Rebhun (University of California, Davis, California, USA). These canine osteosarcoma cells were selected because of their radioresistance (SF2: 0.85) [28], yet they respond well to hyperthermia as a radiosensitizer ( Gy, Gy, and =0.72 Gy with hyperthermia enhancement-factors (EF) and [29] for hyperthermia performed as indicated below.)
Cells were kept in DMEM at 37°C in a humidified incubator with 5% (MCO-18AC-PE, Sanyo, Osaka, Japan). In case of a hyperthermia treatment preceding irradiation, the cells were transferred to another incubator of the same type, set to 42°C, and exposed to a heat-up phase of approx. 40 min, followed by another 60 min of treatment time at the target temperature. The sequence of treatments (hyperthermia followed by irradiation) and the treatment time were chosen to match the clinical practice [30]. To ensure repeatability and quantify thermodynamic effects such as heat transfer and evaporative cooling, incubators were carefully calibrated [29]. In case of an experiment without hyperthermia treatment, the cells remained in the 37°C incubator. Upon completion of the hyperthermia treatment time, the cells were removed from the incubators and irradiated with a 6 MV linear accelerator (Clinac iX, Varian, Palo Alto, USA). Adequate dose build-up and optimal homogeneity of the dose distribution over the irradiation field were ensured by appropriate layers of Plexiglass. Since the irradiation device is also used for regular animal patient treatments, the dose calibration is carried out by a board-certified, qualified medical physicist and is regularly checked with an ionization chamber calibrated at the Swiss Federal Institute of Metrology (METAS).
For logistic reasons (transfer time, setup time, and sequence of irradiation), there was a time-gap of approx. 10 min between the end of the hyperthermia treatment and the beginning of irradiation. Irradiation occurred at doses between 0 Gy and 6 Gy with a dose rate of 600 MU, corresponding to approx. 6 Gy/min. Figure 2 illustrates the timeline of the experiments. It is important to note that while the timeline may suggest otherwise, any experimental procedure (clonogenic and comet assay) discussed below is destructive to the cells. Cells used for a given readout can therefore not be used again for a later or different readout. Thus, the readout originates from different batches of cells.
2.2. Clonogenic Assay
Clonogenic assay is a method to quantify the fraction of cells that survive a treatment, in this case an irradiation event [15]. It works by seeding a number of cells in a dish such that colonies form around these cells due to cell division. After 10 days, the number of colonies are counted and related to the number of cells seeded. If a cell loses clonogenicity due to the treatment, it will not form a clone, while cells which survive the treatment (in the sense of maintaining clonogenicity) will form a colony. The dataset to model clonogenic cell survival in canine osteosarcoma Abrams cells used here was previously published, and the details of the experimental protocol are described in [29].
2.3. Comet Assay
Comet assay is a method to quantify physical DNA damage in individual cells [14] and was performed as follows: approximately 1.5 Abrams cells were seeded in each well of 6-well plates the day before treatments. Cells were treated with radiation and/or heat and harvested after treatments. For this, trypsin was used, and cells were then resuspended in ice-cold PBS. After centrifugation, cells were counted in each sample and resuspended in their DMEM culture medium complemented with 10% DMSO, in an appropriate volume to reach the concentration of 2 cells per mL. Samples for cells used in comet assay were then stored at -80°C. Experiments were repeated 3 times.
Cells from every repeated experiment were thawed on the same day and run for comet assay (5 different runs were needed to run all the samples). After thawing and centrifugation, DMSO was quickly removed and ice-cold PBS added. Cells were suspended in molten LMAgarose (CometAssay® LMAgarose, Trevigen) at a ratio of 1/10 (approximately 1500 cells per sample). Cells were embedded in agarose on a glass slide and left in the dark for 10 min at 4°C. Slides were then immersed in a 4°C lysis solution (CometAssay® Lysis Solution, Trevigen) for 1.5 h in a room at 4°C. Slides were then immersed in the electrophoresis running buffer (8 mg/mL NaOH, 2 mL/mL 0.5 M EDTA pH 8, in O) for 10 min at 4°C in the dark. For electrophoresis, slides were placed in the Trevigen Comet assay tank (CometAssay® Electrophoresis System II, Trevigen) in a cold room, in an exact volume of 850 mL of 4°C electrophoresis solution. Runs lasted 30 min at 21 V and 0.3 A. Care was taken to maintain the same temperature and volume of solution between runs to avoid interrun variability. Slides were finally immersed twice in dHO for 10 min each, then in 70% ethanol for 15 min at room temperature. For staining, diluted SYBR Gold (1 : 10 000, SYBR Gold Nucleic Acid Gel Stain, Invitrogen) was then added to each spot of dried agarose including cells, for 15 min at room temperature, in the dark. Slides were rinsed, dried, and stored at room temperature in the dark.
In order to quantify DNA damage, the microscopy image of the stained comets is analysed with the image processing software COMET IV, which computes a value for each cell/comet, indicating the degree of DNA damage. From the damage metrics offered by the software, the relative tail intensity (RTI) was chosen because it provides a linear relationship between the number of DNA strand breaks and the quantified damage [31, 32]—a property highly desired for the data-mapping introduced in Section 2.5. The resulting data was used in a previous publication [27].
2.4. The Multi-Hit-Repair Model
As mentioned in [27], the MHR model [12] is a dynamic population model where cells are assigned to populations depending on the number of radiation-induced hits (thus the variable name ) they have accumulated. The variable counts the number of cells in population . A hit is defined in this work as a lesion that is hindering the cell from mitosis. In consequence, cells with one or more hits cannot undergo mitosis until all the hits are cured by the repair process. Clonogenicity is the ability of cells to form clones, for which mitosis is a prerequisite. Thus, only the cells in are clonogenic. Figure 3 provides a graphical illustration of the model. Cells can accumulate up to hits, corresponding to the length of the aforementioned chain. The chain length could be infinite, but for an implementation, has to be limited. The practical limit for is chosen such that no congestion at the end of the chain occurs. This criterion is met at ; thus, the chain length was chosen accordingly.
During a simulation run, all cells are clonogenic at first; thus, they are assigned to population , counted by the state variable . Hits are induced by radiation with dose rate , which is set to 0 Gy/min at any time except during irradiation. Thus, is a square pulse that starts at with an intensity of 6 Gy/min (see Section 2.1); the width of the pulse corresponds to the administered dose. While , cells conceptually travel into the chain as they accumulate hits according to a radiosensitivity parameter ( in the context of the MHR model is unrelated to as used in the linear-quadratic model mentioned in Introduction). After irradiation, repair processes inside the cells cure the lesions and thus, cells travel in the opposite direction where they eventually may reach . This repair is governed by the repair rate constant and modulated by a repair function (see below). Alternatively, the repair processes may fail, leading to the death of a cell. This elimination process occurs at a rate of . Thus, the differential equations for population is
DNA repair cannot occur immediately after repair, since radiation not only induces DNA damage but also damages the proteins required for repair. The consequent initial impedance of repair is modeled using the transient biological dose equivalent (TBDE) :
decays after irradiation and is used in the repair function to impede repair after irradiation:
Since some small amount of damage is already present before irradiation, initial conditions were chosen to reflect the damage distribution according to Equation (8) in prior work [27]. Alternatively, it can be assumed that no damage is present before irradiation; i.e., , , and . Negligible differences in terms of the model output were found between these two approaches; thus, the latter, simpler approach is used in this work. The full set of equations is given in the supplementary materials; a summary of the model parameters is presented in Table 1. See [12, 26, 27] for the derivation, validation, and further discussion of the MHR model.
For hyperthermia, the two variables and are introduced to track the state of active () and inactive () repair proteins. These variables represent the respective relative amount of repair protein, and thus, they sum up to 1; i.e., . The activation and inactivation is governed by the following differential equations:
The rate at which inactive repair protein is reactivated, , is assumed to be constant in prior research [12, 26, 27, 30] and throughout this work. While the reactivation may be temperature-dependent, the authors are unaware of any research supporting that hypothesis, thus, following the principle of assuming simple circumstances whenever possible and, is not a function of temperature.
The inactivation of repair protein, however, is temperature-dependent [2, 33, 34]. The rate at which this occurs, , incorporates the Arrhenius law [35] as follows:
The parameter is introduced for numeric reasons. JKmol is the gas constant, and kJmol is the activation energy as published in the literature [12, 36]. It is important to note that may be cell-line specific; thus, the choice of should be revisited in the future once such data is available for the Abrams cell line used here. It is easy to show that the equilibrium of Equations (4) and (5) are and , respectively, for C. Those values therefore serve as initial conditions as it is assumed that this equilibrium is reached prior to the hyperthermia treatment.
The repair function is extended to modulate the DNA repair rate with the amount of inactive repair proteins:
This entails that the rate of repair is reduced both in the presence of inactive repair protein due to thermal effects () and after irradiation when the TBDE is high ().
A temperature of C is set during the hyperthermia treatment. Before and after the treatment, the temperature is set to C.
2.5. Model/Readout Mapping
Since the MHR model is describing radiobiological processes instead of assay readouts, methods need to be implemented to map the model to such readouts. For clonogenic assay, this is relatively straight forward and was introduced in [12]: is tracking clonogenic cells by definition; thus, the surviving number of cells is readily available in . Survival is therefore found by evaluating at the end of the simulation, provided the simulation time is chosen such that the repair process has completed at the end of the simulation.
The mapping to comet data is somewhat more elaborate and was introduced in [27]: The comet readout at a given point in time consists of the quantification of DNA damage in a number of (typically ) cells. Depending on the amount of DNA damage, each cell is assigned to a bin : the first bin tracks the cells with little to no damage and the second bin tracks cells with more damage, etc. The cell count in each bin and population is normalized such that
Finally, the relative bins can be mapped directly to of the MHR model.
In Section 2.4, a hit was defined as an impact on the cell that bars it from mitosis until cured. The correct mapping between physical DNA damage as reported by the comet assay and the model populations presumes knowledge about how much physical DNA damage constitutes one hit. In other words, the relative tail intensities quantifying DNA damage must be scaled prior to the mapping to to maintain the semantics implied by the MHR model (i.e., the definition of a hit). In [27], the correct scaling factor was unknown, and thus, tail intensities between 0 and 4% were mapped to arbitrarily (as discussed there, the model can still be used with a wrong scaling factor, but the parameters may lose the meaning they were originally designed for). In this work, the scaling factor is not fixed to a single, convenient value arbitrarily. Instead, the procedure is repeated with different scaling factors within a sensible range. In order to achieve this, the scaling is formalized by the variable which denotes the largest tail intensity that is still mapped to . Hence, in the above example.
Interestingly, the method failed to reproduce experimental comet data for small values of . For large values of , the resulting parameter values were in violation of the lower bound stipulated by Equation (9) (see Figure S1). Only a small region around was free from these issues; thus, was used.
2.6. Approximate Bayesian Computation
Approximate Bayesian computation (ABC) [37] is used to estimate distributions of model parameters. The method works as follows: for each parameter, the range of biologically meaningful parameter values is estimated. For example, with h (the upper boundary of this parameter), repair proteins reactivate very quickly from the irradiation event; the repair probability recovers to 94% 30 min post irradiation. This is unrealistically high given the typical delays observed experimentally (see Figure 4 and [38]). Since no prior information is available on a given parameter values’ positions within the search space , uniform prior distributions with boundaries and , , are chosen. The boundaries are listed in Table 1 for each parameter. Determining the lower bound on presents a special case: it is easy to show that in the absence of any repair (i.e., ),
for the duration of irradiation. After irradiation, , and thus, remains constant. Because is mapped to the fraction of surviving cells (see Section 2.5), a lower bound for can be established by solving Equation (9) after substituting for as reported in the clonogenic assay and setting to the point in time at which irradiation ends.
At the beginning of the parameter search, sets of parameters are initialized by drawing from the prior distributions. Predictions are made by running the model in a forward fashion, extracting the predicted readout as described in Section 2.5 and comparing it to experimental data. This yields an error according to Equations (10) and (11) (see below).
In each iteration of the search, the parameters are perturbed; the new parameter values are kept if decreases and are discarded otherwise. In a simulated annealing fashion [39], the amount of perturbation is gradually decreased as the search progresses. A cut-off value of was chosen. In the end, sets of parameters are left; all of which provide a satisfactory error. In this work, was chosen with 250 iterations.
The objective function for the parameter search with cell survival data is
for the radiation doses , the experimentally obtained surviving fraction of cells , and the predicted surviving fraction of cells .
Similarly, the objective function for the parameter search with comet data is
for time point , normalized population , and normalized comet readout as defined in Equation (8). A combined calibration was attempted with a combined objective function (see discussion in Section 4).
2.7. Software
The methods discussed above were implemented in python (version 3.5.2) using the abcpy module [40] for ABC (version 0.5.3). R version 3.6.0 was used to create the plots; the code and data are available online (https://github.engineering.zhaw.ch/weyl/synthetic_comet). The software can be configured to use either input data from clonogenic assay or input data from comet assay. Depending on this selection, the corresponding objective function or is used. The results shown in Figure 5 are obtained with the software running on clonogenic mode, i.e., evaluating , while those in Figure 4 are obtained with the software running in comet mode, i.e., evaluating .
3. Results
Two model outputs for cell survival are shown at the top of Figure 5, as produced by the software running in clonogenic mode. The examples were chosen according to similar : for both instances, . Experimentally, it would be very challenging (if not impossible) to discriminate between the two curves. Yet, the parameters and the dynamics shown at the bottom are very different from each other: in the left case, most hits have vanished after 2 h, while the same requires 4 h in the right case.
With the software running in comet mode (i.e., minimizing ), results are shown in Figure 4. The data on the left represents a random pick from the parameter result set and produces a cell survival curve very different from cell survival found experimentally (=0.12). The ones on the right is the curve with the lowest error found in the set ().
Figure 6 shows a histogram panel of the parameters unrelated to hyperthermia (see Figure S2 for parameters , , and ). In the top row, parameters from the software in clonogenic mode are shown while in the middle row, parameters from the software in comet mode are shown. The bottom row shows the joint distribution, calculated from the previous two rows. Generally, values for and are centered around one or two peaks, while, e.g., is more uniformly distributed in the comet case, but clonogenic assay data suggests that the parameter peaks at low values.
In addition to joining the two posterior distributions for each parameter, a calibration was attempted where the two objective functions were combined with a weighting factor :
In order not to prefer any assay source from the other, errors from the previous single-assay runs were used to select such that the two terms are of the same order of magnitude. This attempt failed; the ABC solver never left its seeding state (see discussion in Section 4).
4. Discussion
The results shown in the previous section clearly call for a combined approach, where both clonogenic assay and comet assay data are used as sources of information for parameter search. However, the traditional approach of combining two objective functions failed. This is because in the seeding state, ABC with simulated annealing rejects samples from the prior that are above a certain threshold (values up to were tried). Since it is difficult to find parameters that satisfy both objective functions, the seeding state never completed. Thus, the computationally much lighter approach with joint posteriors is proposed, allowing for additional flexibility in combining further calibration results.
The results in Figure 5 reveal that survival curves lack sufficient information for MHR model calibration. This is hardly surprising, as it was argued before that the clonogenic assay captures information very distant from the process that is being modeled. On a side-note, any attempt to calibrate a model with 8 parameters from 5 data points is likely going to fail, which is yet another reason to include additional data sources. However, even with this little information, the top row in Figure 6 reveals regions of interest for some parameters, e.g., for and . On a related note, the resulting posterior distributions for and are bimodal. This is an important finding that is concealed by a method aiming at point-estimates, such as differential evolution one used in [27]. Indeed for the value, one peak of the histogram corresponds to the range of parameters found in [12], while the other peak corresponds to the range of parameters found in [27]. Interestingly, these two regimes also correspond to the two instances depicted in Figure 5.
Based on the aforementioned rationale, one may assume that the use of comet assay readouts would cure these issues. However, Figure 4 demonstrates that this is not the case. Otherwise, any parameter set would yield an adequate cell survival curve. Discussing potential explanations for this observation is critical since the resulting conclusions govern the choice of further data to address the open issues: for quantification of the damage, the relative tail intensity is assessed for approx. 100 cells per assay. This quantification does not discriminate between cells that have a chance to reach , cells that have already initiated apoptosis and will never reach , and cells that are on the brink of death for other reasons. In fact, the quantification may even contain cells that are already dead but still have DNA that is visible in the microscopy image. However, the ability to reach is essential for producing a survival curve from the model. Thus, a possible explanation for the inability to achieve successful model calibration from comet assay readout alone could be that the readout does not carry sufficient information about the ability to reach . Furthermore, dead cells that have degraded so far as to not have any quantifiable DNA whatsoever would not be considered for comet assay, and the normalization of the 100 cells to a relative histogram would mostly masquerade their existence: the only way for dead cells to influence the results is in the ratio , since a surviving cell would contribute to (thus increasing ), but if the same cell had died, it would not contribute to any (thus increasing for ).
For the parameters , , and (Figure 6) as well as the parameters related to hyperthermia (Figure S2), uniform posterior distributions are obtained. The method thus reveals that more input data is required to identify these parameters. Parameters and relate to a transient repair inability due to the irradiation event. From Figure 4, it can be seen that this effect vanishes approx. 30 min. after irradiation. Thus, further data within that time frame could yield better estimates for those parameters. could be identified by running a series of clonogenic assays at various dose rates. At low dose rates, irradiation would not be considered as an event but have a finite duration and repair may start already during irradiation. Such dose-rate-dependency was shown in the past to be reproduced by the MHR model [12], and the rate of repair could become identifiable. The hyperthermia parameters could be refined with data from a study with varying time-gaps. Such data from clonogenic assay has been publishedI [9] but not from comet assay and with different cell lines. As mentioned in Introduction, assessing the repair-protein reactivation rate constant would be of great clinical use, as it would allow a better assessment about tolerable time-gaps (and variation thereof) between irradiation and hyperthermia. Since the order of the two treatments (i.e., hyperthermia prior to versus after irradiation) was shown to have minimal effect on cell survival [9], additional input data in this regard would likely not improve the calibration results. Investigating more cell lines would reveal which parameters may vary by how much between subjects.
Clonogenic cell survival and comet assay measurements were shown to be repeatable [14, 29]; thus, it is reasonable to expect repeatable results from patient biopsies [38]. This would allow for a per-patient calibration, e.g., to improve the treatment plan on a per-patient basis. As mentioned in the previous paragraph, such an endeavour would require data from appropriate sources to identify the relevant parameters, rather than just more amounts of data.
In case the model cannot be calibrated at all despite these efforts, it could be simplified, for example by replacing the TBDE with a fixed window of no repair after irradiation, removing the parameters and . Alternatively, it is conceivable to split the process at the time of the irradiation, yielding a hyperthermia process that sets up the initial conditions for a subsequent DNA damage and repair process. Splitting the model in this way could yield closed-form solutions or approximations thereof for some state variables, paving the way for a much simpler calibration strategy.
The model and the strategy presented in this work have a number of potential limitations. First, the model does not incorporate any mitosis, which occurs without doubt in until the cells are fixed and the clonogenic assay is performed. However, one can argue that for a given cell line, any mitosis would occur at a fixed rate. While the number of cells would increase, their ratio would remain the same. Because the clonogenic cell survival assays used in this work are normalized, mitosis cancels out. Some cells may, however die only after a few cell cycles. This falls in the gap between the last comet assay and the point in time when clonogenic assay is performed and is not modeled in the MHR model. Second, the model does not incorporate any effects of direct cytotoxicity, i.e. thermal cell-killing. This is alleviated by the fact that such direct cytotoxicity was not observed in any of the control experiments performed with hyperthermia alone [29]. Third, the model does not correctly describe inhibition of DNA repair proteins above a temperature threshold of 42.5°C–43°C, since those proteins are believed to enter a different regime above that threshold [35]. While this is a limitation, it does not affect the work presented here since the highest temperature applied in vivo and in silico was 42°C.
5. Conclusions
This work demonstrates that a holistic approach is necessary to calibrate the MHR model parameters. Relying on clonogenic assay data or comet assay data alone, as it has been done in the past, proved to be insufficient to establish unambiguous model parameters. Even with this combined approach, some parameters remain unidentified. However, the ABC method has the advantage of joining existing posterior distributions with distributions obtained from calibration runs with new input data. This ability is critical since model calibration with ABC is, despite all its advantages, very slow. Combining posterior distributions from ABC, however, is fast. Following this approach, data from different assays can be combined in a modular fashion without the need of rerunning the full calibration.
While the application of the method presented is radiotherapy, hyperthermia, and treatment planning, the method presented here addresses a more general problem; thus, many other instances exist where the application of this method would be of value.
Data Availability
The data are available online (https://github.engineering.zhaw.ch/weyl/synthetic_comet) along with the software code required to perform the analyses shown in this work.
Conflicts of Interest
The authors declare that there are no conflicts of interest.
Acknowledgments
This work was supported by the Swiss National Foundation (grant number 320030_163435)—Stephan Scheidegger, Carla Rohrer Bley. This work was financially supported by a grant from the Marie-Louise von Muralt-Stiftung für Kleintiere.
Supplementary Materials
The supplementary materials contain (A) the full equation system of the multi-hit-repair (MHR) model and (B) supplementary figures with additional parameter histograms. (Supplementary Materials)