A cellular automaton model for spheroid response to radiation and hyperthermia treatments

Thermo-radiosensitisation is a promising approach for treatment of radio-resistant tumours such as those containing hypoxic subregions. Response prediction and treatment planning should account for tumour response heterogeneity, e.g. due to microenvironmental factors, and quantification of the biological effects induced. 3D tumour spheroids provide a physiological in vitro model of tumour response and a systems oncology framework for simulating spheroid response to radiation and hyperthermia is presented. Using a cellular automaton model, 3D oxygen diffusion, delivery of radiation and/or hyperthermia were simulated for many (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bf{1}}{{\bf{0}}}^{{\bf{4}}}{\boldsymbol{-}}{\bf{1}}{{\bf{0}}}^{{\bf{7}}}$$\end{document}104−107) individual cells forming a spheroid. The iterative oxygen diffusion model was compared to an analytical oxygenation model and simulations were calibrated and validated against experimental data for irradiated (0–10 Gy) and/or heated (0–240 CEM43) HCT116 spheroids. Despite comparable clonogenic survival, spheroid growth differed significantly following radiation or hyperthermia. This dynamic response was described well by the simulation (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\bf{R}}}^{{\bf{2}}}$$\end{document}R2 > 0.85). Heat-induced cell death was implemented as a fast, proliferation-independent process, allowing reoxygenation and repopulation, whereas radiation was modelled as proliferation-dependent mitotic catastrophe. This framework stands out both through its experimental validation and its novel ability to predict spheroid response to multimodality treatment. It provides a good description of response where biological dose-weighting based on clonogenic survival alone was insufficient.

A cellular automaton model for spheroid response to radiation and hyperthermia treatments Sarah c. Brüningk * , peter Ziegenhein , ian Rivens, Uwe Oelfke & Gail ter Haar thermo-radiosensitisation is a promising approach for treatment of radio-resistant tumours such as those containing hypoxic subregions. Response prediction and treatment planning should account for tumour response heterogeneity, e.g. due to microenvironmental factors, and quantification of the biological effects induced. 3D tumour spheroids provide a physiological in vitro model of tumour response and a systems oncology framework for simulating spheroid response to radiation and hyperthermia is presented. Using a cellular automaton model, 3D oxygen diffusion, delivery of radiation and/or hyperthermia were simulated for many (10 −10 4 7 ) individual cells forming a spheroid. the iterative oxygen diffusion model was compared to an analytical oxygenation model and simulations were calibrated and validated against experimental data for irradiated (0-10 Gy) and/or heated (0-240 CEM 43 ) HCT116 spheroids. Despite comparable clonogenic survival, spheroid growth differed significantly following radiation or hyperthermia. this dynamic response was described well by the simulation (R 2 > 0. 8 5 ). Heat-induced cell death was implemented as a fast, proliferation-independent process, allowing reoxygenation and repopulation, whereas radiation was modelled as proliferation-dependent mitotic catastrophe. this framework stands out both through its experimental validation and its novel ability to predict spheroid response to multimodality treatment. it provides a good description of response where biological dose-weighting based on clonogenic survival alone was insufficient.
The combination of radiotherapy (RT) with hyperthermia (HT, prolonged elevation of temperature above the physiological range) provides a promising treatment for radiation-resistant tumours by invoking a mechanism for localized heat-induced radiosensitization, in cases for which ionising dose escalation is limited due to normal tissue toxicity. The biological effects induced by RT and HT act synergistically [1][2][3][4][5] , but are strongly dependent on the heating temperature, duration and sequence of the combined treatment. Quantification of the biological effects produced is essential for treatment dose prescription and response prediction. Biological response quantification based on clonogenic cell survival data has been proposed [5][6][7][8] with dose weighted by a relative biological effectiveness using heating time and temperature dependent parameters to calculate the biologically equivalent dose (BEQD) in units of RT only treatments. This type of biological effect calculation does not, however, account for the dynamics of tumour growth delay or shrinkage following treatment which is affected by the treatment modality due to differences in the cell death mechanisms induced 9,10 . In addition, it has been shown that cells cultured in 2D, as used for clonogenic assays, may respond differently to treatment than more physiologically relevant 3D tumour spheroid cultures [11][12][13] raising concerns about the applicability of biological dose weighting purely based on clonogenic cell survival assays. Tumour spheroids are avascular 3D aggregates of cells which mimic the physiological microenvironment of tumours in vivo more closely than 2D cultures by providing gradients of oxygen, nutrients and waste products. Cells are in direct (3D) contact with each other, thus stimulating further inter-cellular communication and allowing interaction between cells and their extra-cellular matrix. Tumour spheroids are therefore considered to provide a simple in vitro tumour model that allows study of treatment-induced growth dynamics of the cell population as a whole, and may be a better in vitro model for quantifying and predicting in vivo tumour growth than 2D cultures [14][15][16] . Recent studies have investigated the dynamics of cell cycle progression and cell death within irradiated spheroids 9,10 and suggested a proliferation dependence of radiation-induced cell death as a consequence of irreparable DNA damage and mitotic catastrophe 17,18 . Although investigations of spheroid response to HT or combined RTHT treatments have been performed [19][20][21][22][23][24] , these mainly focused on an analysis of spheroid growth delay, cell viability and dissociated spheroids clonogenic survival. An increased heat resistance of cells treated in spheroids has been suggested. Little information exists, on the dynamic process of heat-induced cell death and to what extent differences in mircoenvironmental factors and cell cycle distribution contribute to the observed increase in heat resistance. We suggest that within heated spheroids cells undergo a rapid, and proliferation-independent cell death. This agrees with the observation that several cellular components, including functional and structural proteins, are affected by this treatment 25,26 . Current BEQD models cannot account for these differences in the cell death mechanisms leading to microenvironmental changes and therefore are unable to explain the observed differences in spheroid growth following treatment with RT or HT at isoeffective (thermal) dose levels. Agent-based systems oncology simulation, such as cellular automaton models [27][28][29] , provide the mathematical tools required to meet this objective. Here, cells are modelled as individual voxels on a (3D) computational grid with specific growth and treatment response properties. Dynamic changes in microenvironmental factors (e.g. oxygenation) are also accounted for. We present a 3D cellular automaton (CA) model for simulating simultaneous RT and HT treatments, and compare it to experimental data obtained using spheroids formed from the human colorectal carcinoma cell line HCT116. The aim of this study is to demonstrate the importance of considering differences in the cell death mechanisms induced after RT or HT treatments and how these impact the dynamic growth response of the spheroid as a first step towards more sophisticated simulations and BEQD estimations for tumours in vivo. Where possible, a deliberately simple implementation was chosen in order to minimize the number of model parameters.

Methods experimental data.
Simulations were compared to experimental data obtained for the colorectal carcinoma cell line HCT116 that will be published separately and describe in detail spheroid initiation, assessment of their growth, live cell imaging and preparation of histological sections from spheroids, as well as delivery of radiation and hyperthermia treatments for spheroid and clonogenic survival analysis. Briefly, cells were grown in Dulbecco's modified Eagle's medium (Gibco Life Technologies Ltd. Paisley, UK) supplemented with 10% fetal bovine serum (PAN Biotech GmbH, Aidenbach, Germany) under standard mammalian cell culture conditions (37  C, 5% CO 2 in a humidified atmosphere). Spheroids were prepared by seeding 300 HCT116 cells in exponential growth phase, in U-shaped, ultra low attachment plates (ULA plates, #7007, Corning B.V. Life Sciences, Amsterdam, NL), and were incubated for four days prior to treatment. This allowed the cells in each well to form a single dense 3D aggregate, and gradients of nutrients and oxygen to be established. Spheroids were irradiated in a small animal radiation research platform (220 kVp, 13 mA, 0.14 mm Cu filtering, X-Strahl, Camberley, UK) and/or heated in a Biorad Tetrad2 DNA Engine PCR thermal cycler (Bio-Rad Laboratories Inc., Berkeley, USA) at 47  C for various lengths of times, before being transferred in complete growth medium supplemented with 1% penicillin, streptomycin and amphotericin B (all Sigma Aldrich Ltd., Dorset, UK) and 5 µg ml / propidium iodine (PI, Sigma Aldrich -which stains dead cells) to fresh ULA plates. Cultures were maintained for three weeks after treatment with the medium being renewed every three days. Spheroid diameter (growth) was assessed every other day using an automated Celigo TM imaging cytometer (Nexcelom Biosciences LLC, Lawrence, USA). An implementation of a deep learning neural network was used for automated spheroid segmentation. computational model. A previously published cellular automaton framework 30 was extended to describe spheroid growth on a 3D lattice. A typical grid size was 200 × 200 × 200 cubic voxels, modelled with 0.012 mm edge length, corresponding to the average size of HCT116 cells obtained previously for 2D simulations 30 . Cells are initialized as forming a dense sphere on the grid, with each cell occupying one voxel. Cells progress through a predefined cycle (G1, S, G2, M-phases) according to an individualized timer and starting from a random cycle stage as described previously 30 . During cell division, next neighbours were searched alternately in the 3D Moore and vonNeumann neighbourhoods up to second order. Here, more central locations are favoured to ensure a compact, spherical structure. If no free space was available, cells entered a reversible quiescent stage (G0) until neighbouring spaces are vacated. In the absence of treatment, cells were simulated as having an infinite lifespan and unlimited division potential. However, oxygenation levels and the presence of free next neighbours influenced cell cycle progression and cell survival (death after prolonged hypoxia) as described in the following section. No cell migration within the spheroid was considered.
Oxygenation and untreated spheroid growth. The partial oxygen pressure p t x ( , ) O2 at 3D location x and time t were modelled by solving the oxygen reaction diffusion equation iteratively: Here, D O 2 is the oxygen diffusion coefficient in the medium, and Φ x ( ) denotes the cellular oxygen consumption rate where voxel x is occupied (Φ = x ( ) 0 for unoccupied voxels). Boundary conditions are set to full oxygenation outside the spheroid volume corresponding to a partial oxygen pressure of 100 mmHg as previously reported by Grimes et al. for spheroid cultures within a CO 2 incubator 31,32 . Since the time scale for oxygen diffusion (seconds) is far less than that of cell cycle progression (hours), the oxygen distribution is updated within each cell cycle iteration until an equilibrium stage is reached. Equilibrium was defined as a maximum change in oxygen partial pressure of ≤0.01% at any location.
It is known that cells die if kept under hypoxic conditions for prolonged periods of time (chronic hypoxia) 33 . Lack of other nutrients, accumulation of cellular debris, and mechanical pressure are known to contribute to cell cycle arrest and cell death especially at the core of a spheroid [34][35][36][37] . For simplicity, only hypoxia-induced cell death was modelled in this framework, but the parameters may be regarded as accounting for several microenvironmental factors summarized under the term "hypoxia".
www.nature.com/scientificreports www.nature.com/scientificreports/ Cells were labelled as hypoxic once the partial oxygen pressure fell below 11 mmHg 38 . These hypoxic cells died with probability p hypoxiaDeath at each simulation time step if they did not return to normoxic conditions, thus forming a necrotic core. Although the underlying mechanism is still subject to current research, it is agreed that the formation of a necrotic core results in spheroid growth plateauing over time. To account for this experimental observation in the simulation, necrotic cells were removed from the computational grid with probability p clearNecrotic at each time step, whereas a constant growth of the proliferating spheroid layer was assumed. All cells were continuously shuffled inwards in order to maintain a dense spherical structure. By choosing small values for p hypoxiaDeath and p clearNecrotic it was not necessary to define additional parameters to control the minimum time under hypoxia before cell death, or a minimum time after cell death before clearance. This implementation clearly does not reflect the actual biological processes, instead it provides a simplification that mimics the spheroid growth plateau while being restricted to a fixed-sized lattice.
Radiation response. Radiation-induced cell death in 3D spheroids was modelled as a proliferation dependent process. Cells were randomly divided into a dying, and a surviving population according to the clonogenic surviving fraction S RT that was calculated using the linear quadratic model 39 with cell line specific parameters α and β as described previously 7 . Weighting factors for cell cycle stage, γ, and oxygenation status, d OER , were applied upon calculation of this surviving fraction.
Here, the decreased radiosensitivity of hypoxic cells relative to normoxic ones was accounted for by weighting the radiation dose, d, by the oxygen enhancement ratio, OER, to obtain the apparent delivered dose, d OER .
As the ioseffective ratio of doses under normoxic and hypoxic condion, the oxygen enhancement ratio (OER) 40 depends on the partial oxygen pressure pO 2 . However, the specific cell line, and linear energy transfer of the radiation also contribute to the OER. Typical experimentally determined values for the OER are in the range 1-3 32 . This motivated the choice of low and high oxygenation level limits of the step function described in equation 3 using a linear interpolation for hypoxic oxygen levels. Figure 1A,B shows a schematic diagram and the probability driven response cascade used to mimic radiation-induced cell death. Cells randomly selected to survive the radiation treatment (with probability S RT ) continued to proliferate, no cell cycle delay was included for surviving cells. "Dying cells" underwent mitotic cell death with probability p mitoticDeath each time they attempted division at the end of M-phase; with probability − p 1 mitoticDeath that an irradiated cell successfully divided into two daughter cells, both labelled as dying. Experimentally, irradiated spheroids showed no signs of radiation-induced growth inhibition within the first days after exposure. This was accounted for in the simulation as a two stage process using two parameter values for p mitoticDeath : within the initial days (parameter t delayRT ) < . p mitoticDeath 0 5 , i.e. net growth, and > . p mitoticDeath 0 5 , for the remaining time, resulting in spheroid shrinkage.
The implementation presented allows gradual shrinkage of irradiated spheroids from the outside inwards. Quiescent cells were sequentially reactivated and underwent mitotic death once attempting division. This resulted in a continuous shrinkage of a lethally irradiated spheroid (i.e no surviving cells), while maintaining a compact spherical structure. If the proportion of surviving cells was sufficiently large, spheroid regrowth took place before the inner-most quiescent cells re-enter the cell cycle. These dying quiescent cells remained within the spheroid core, encapsulated by layers of proliferating, surviving cells. Figure 1A illustrates the dynamic cell death after the irradiation described above for a high (very few surviving cells) and low (few dying cells) dose.
Hyperthermia treatments. For simulation, heated cells were separated into surviving and dying populations according to the clonogenic surviving fraction, S HT , directly after heat treatment. Since the thermal-resistance observed in HCT116 spheroids was enhanced compared to their clonogenic survival in 2D monolayer cultures, adjustments were made to the clonogenic surviving fraction calculated using the AlphaR model 7 . It was assumed here that the 3D surviving fraction in response to heat treatment, S t ( ) HT 43 , followed a simple step-wise function of thermal dose, t 43 , the 2D culture AlphaR model parameters β HT and α 0 , and the plateau surviving fraction S HT plateau HT HT HT plateau 43 43 43 43 www.nature.com/scientificreports www.nature.com/scientificreports/ No adjustments were made to allow for differences in cell death as a function of the cell's oxygenation level as previously suggested 41,42 . Although other parameters, such as a more acidic pH in the spheroid core, would enhance heat-induced killing, these were not considered in this simplified implementation. Cell cycle stage dependent heat sensitivity variation was considered, in particular differences between cycling and quiescent (G0) cells. This was achieved by using a weighting factor (γ HT ) as described above for the modelling of radiation response. A constant ratio of 1.5 (as used in previous 2D simulations 30 ) between G1, S, G2-phase cells was used. G0 cells were assumed to be more resistant by a factor of three relative to the average surviving fraction of the cycling cells, motivated by 43,p.19 .
Heat-induced cell death dynamics were implemented as occurring after a delay time (sampled from a normal distribution, mean value t delayHTtodeath ), followed by immediate clearance of dead cells. Surviving cells were simulated as continuing their cell cycle progression after a randomly assigned cycle arrest (t cellcycledelay ), sampled from a uniform distribution. Free parameters for the HT response cascade were therefore S HT plateau , , t delayHTtodeath , and t cellcyclearrest . Figure 1D shows the response cascade implemented to model the behaviour described.
Since cells are removed from the computational lattice independent of their position, this model allows reoxygenation and re-activation of previously quiescent cells. Due to the nature of the culture wells used (U-bottom shape), cells are, however, drawn towards the centre of mass resulting in a re-shuffling of the surrounding cells to fill any gaps. Figure 1C shows schematic images of heat-induced cell death mechanisms which resulted in a destabilization and disintegration of the spheroid for large thermal doses (few surviving cells), before the spheroid re-formed from the surviving cells. Parameter fitting and validation procedure. In the CA model, the diameter of the simulated spheroid was calculated as the average of the 100 maximum distances of cells to the spheroid centre in MATLAB (MathWorks, version 2017a). The coefficient of determination, R 2 , was used as a measure of agreement between simulated and experimentally measured growth curves. Simulation parameters were adjusted to maximize R 2 in the calibration step using a grid search. Since the number of free parameters was limited as much as possible, each parameter controlled a specific aspect of the model and could thus be independently fitted to the specific data using a "one at a time approach". Below, the specific features used for fitting each parameter are given in brackets. Since parameter cross correlations were not applied during the fitting procedure, multiple parameter combinations yielding a similar solution were not encountered in this approach. It should be noted that absolute parameter values refer entirely to an artificial, simplified biological system and thus should not be used for direct interpretation of biological meaning.
Cell growth was modelled with an allowance of 3 days for cell attachment and acclimatization relative to experimental growth curves. I.e. simulation started on day 3 in terms of the experimental time since seeding the cells and all growth curves are shown with the relevant time correction.
Oxygenation and untreated spheroid growth parameter estimation. Oxygen distributions within spheroids were compared qualitatively to histological sections stained with pimonidazole (for hypoxia), DAPI (for cell nuclei to show cell distribution) and Ki-67 (to show proliferation). For quantitative evaluation, simulations were compared with an analytical oxygenation model proposed by Grimes et al. 31 for which parameters had previously been estimated for HCT116 spheroids 44,45 . This model provides a description of the partial oxygen pressure, p O2 , as a function of the radial distance to the spheroid centre, r, based on experimentally measurable dimensions within stained histological spheroid sections. Here, p 0 denotes the partial oxygen pressure in the culture medium ( ≈ p mmHg 100 0 31 ), r l is the diffusion limit, r n is the radius of the anoxic necrotic core, and r 0 is the outer spheroid radius. The diffusion limit r l , i.e. the minimum radius of a spheroid with partial oxygen pressure equal to zero at its centre, can be measured experimentally since it is directly related to the width of the proliferating and quiescent zone, r c , and the spheroid radius r 0 . Using these values, the partial oxygen pressure as a function of the radial distance to the spheroid surface was calculated using equation 6. Results were compared with the simulated oxygen distribution using an oxygen consumption rate of 22.1 mmHg/s as proposed by Grimes et al., and the oxygen diffusion coefficient D O 2 (gradient of hypoxia) was adapted to minimize the difference between central cross-sections of simulated and analytically calculated oxygen distributions.
Having calibrated the oxygen diffusion model, the parameters controlling untreated spheroid growth, i.e. initial number of cells (initial speroid diameter), cell doubling time (growth rate under full oxygenation), p hypoxiaDeath (onset of growth plateau), and p clearNecrotic (slope of growth plateau), were adjusted to maximise the agreement (in terms of coefficients of determination) between simulated and experimental growth curves of untreated HCT116 spheroids.
Treatment response parameter fitting and validation method. The experimental data for 10 Gy were used to determine the fitting parameters by adapting p mitoticDeath , before and after an initial delay (slopes of the curve), t delayRT (peak location of the curve), to cell death after irradiation (free parameters: p mitoticDeath 0 5 , t delayRT ). The 10 Gy growth curve was chosen, since spheroids did not regrow after this treatment within the three week observation period. This simulation was not subject to uncertainties of surviving cells trapping dying G0 cells which simplified parameter fitting. After parameter fitting, growth curves after 2 Gy and 5 Gy irradiation were simulated without further parameter changes (validation). Similarly, free parameters S HT plateau , (time of regrowth), t delayHTtoDeath (rate of shrinkage), and t cellcycle (onset of shrinkage), were fitted for HT treatment based on the 240 CEM 43 growth curve and these parameters were then used to predict 40 and 80 CEM 43 treatments. For combination treatments, no further changes to the simulation parameters were made and the growth curves after treatment with 2 Gy + 40 CEM 43 , 2 Gy + 80 CEM 43 , and 2 Gy + 160 CEM 43 were simulated and compared to experimental data.  www.nature.com/scientificreports www.nature.com/scientificreports/ levels are shown in Fig. 4. Simulated and analytically calculated distributions agreed with differences of less than 15%, the maximum deviations being observed in the outermost cell layers.

Results oxygen distribution and untreated spheroid growth.
Model parameters controlling the formation of a necrotic core leading to a growth plateau were fitted with the best fit ( = . Simulation of treatment response. Parameter fitting to the 10 Gy growth curve was optimal (R Gy thereafter. The results obtained, together with the validation data sets for 2 Gy and 5 Gy irradiated spheroids are shown in Fig. 6 (top). Results are shown for a representative simulation run using the mean surviving fraction (see Eq. 2). For 2 Gy and 5 Gy treatments, simulations were also performed using the mean ± one standard deviation of the clonogenic surviving fraction to show the variations in the data. Since 10 Gy was considered to be a treatment that killed all cells within the spheroid, no variation in the surviving fraction was accounted for. There was good agreement between simulated and experimentally measured data, for both the calibration and validation data sets (R Gy   ), or to a breakdown of the extracellular matrix which was not accounted for in the simulation. It was, however, considered that correct prediction of spheroid growth delay was most important, as was achieved (see Fig. 6 (middle)).
The simulation of combined RTHT treatments was carried out without further modifications of the parameters used. Simulated and experimentally measured spheroid growth curves are shown in Fig. 6 (bottom). There was good agreement for this data set (

Discussion
Modelling spheroid response in terms of volumetric growth required close mimicking of the biological response observed experimentally. Despite using a deliberately simple implementation of a CA approach, it was possible to reproduce the most important features of the treatment of 3D tumour spheroids with radiation, hyperthermia, or both in combination. These included modelling oxygenation to predict the formation of a necrotic core and spheroid growth plateau, as well as implementing proliferation dependent radiation-and proliferation independent heat-induced cell death. It was shown that the simulation modelled experimental spheroid growth response to RT, HT and RTHT treatments (Fig. 6) (all > . R 0 85 2 ) accurately. Our framework represents an improvement over previous similar models 46,47 , both through the extent of experimental validation performed and the different treatment modalities simulated, i.e. RT, HT, and RTHT. Sophisticated biological interpretations 9,48 and computational models 46,47 have been proposed previously for the modelling of radiation response as a single modality treatment or in combination with chemo therapy. We identified and implemented the underlying cell death mechanisms between radiation and heating, and subsequently determined their influence on spheroid growth. We suggest that similar differences in cell death mechanisms are likely present for radiation and drug treatments due to different cellular targets, and it may be essential to consider dynamic cell death for such multi-modality therapies, to allow calculation of BEQD distributions for any kind of multimodality therapy.
One of the aims of the experimental calibration was to facilitate modelling of oxygenation effects. These include the formation of a necrotic core following prolonged hypoxia, the plateauing of spheroid growth, and the adaptation of the oxygen enhancement factors to capture the influence of hypoxia and normoxia in irradiated spheroids. The HCT116 spheroids used here were a good model for adapting the relevant response cascade as spheroids were well oxygenated on treatment day but developed a hypoxic and necrotic core over the observation period (three weeks). Although reference values were available for the oxygen consumption and diffusion rates within HCT116 spheroids 44 , using these values directly would imply that the voxel size used matched experimental cell sizes. Therefore, a small adaptation of the relevant diffusion rate was made, based on comparison with the analytical model proposed by Grimes et al. for modelling the partial pressure of oxygen. Despite this adaptations, differences of up to 15% in partial oxygen pressure were observed between the two approaches, with the analytical calculation providing higher values within the proliferating spheroid zone. Extreme values of 15% may be explained by a some individual cells lying outside the calculated mean diameter of the spheroid, whereas the analytical model assumes a perfect sphere. Although there was only a relatively homogeneous off-set between distributions on day four, larger spheroids showed only small differences with respect to the oxygenation gradients. These may originate from differences in the "growth history" of the spheroids used in the model by Grimes et al. and the ones used here. Grimes et al. used the liquid overlay technique for spheroid production and used spheroids 4 to 17 days after formation for their analysis and parameter estimation. This implies a larger range in incubation times, and a wider range of spheroid sizes obtained than the more standardised time points and spheroid formation method used here. There was qualitative agreement between simulated and experimental oxygenation levels estimated from histological sections, showing little significant hypoxia in small (day 4) spheroids, with increasingly hypoxic regions for larger spheroids (days 11 and 25). Parameters controlling the influence of hypoxia on treatment response, such as the OER, could not be fitted here, due to the absence of significant hypoxia on treatment day (see Fig. 4). Spheroids with hypoxic, but not necrotic, cores should be evaluated for this purpose. This could be achieved by treating the same HCT116 spheroid at a later time point after seeding.
Other than oxygen distribution, the current framework did not include modelling of other diffusive factors, such as glucose, waste products, or pH-levels that have been reported to affect heat sensitivity of monolayer cultures and cell viability in general and are therefore likely to affect treatment response of cells within a spheroid 49,50 . It should therefore be noted that the current formulation of "oxygenation" may indeed emcompass a number of diffusive processes to limit the number of model parameters. Accounting for separate description of oxygen and other diffusive factors could be the subject of future model developments if deemed necessary.
Another limitation of the frameworks was the type of simulation framework chosen, a cellular automaton model which intrinsically limited the ability to model some potentially relevant biological processes. Due to the fixed lattice size, it was not easy to implement treatment-induced variations in cell size, a loosening of the cellular structure (i.e. cell density variation), or mechanical compression of cells within the necrotic spheroid core. A potential solution to this problem would be the use of a lattice-gas model 29,51 , providing multiple cell occupancy for each constant sized voxel, independent of the number of cells present. This would allow modelling of varying cell density and simplify the simulation of cell migration that could not be considered here in addition to the constraint of a dense spherical structure. Despite this limitation, our model was able to describe the observed spheroid growth curves well, and we thus deemed it sufficient for the current purpose.
The ultimate goal for the simulation framework will be modelling of vascularized tumours in vivo and the simulation presented here represents a first step towards this. Using a systems oncology model based on the one presented here, calibrated for each hierarchically relevant environmental detail (i.e. cellular properties, inter-cellular interactions, 3D growth and microenvironmental influences) would allow comparison of the importance of intrinsic cellular properties, such as growth rate, cell cycle progression, and clonogenic treatment response, with the influence of microenvironmental factors such as hypoxia, pH, etc. This should help bridge the gap between experimental results observed in vitro in 2D or 3D cultures and those in vivo. It is generally agreed that hypoxic cells are more radioresistant than normoxic ones when evaluated by clonogenic assays. In an in vivo tumour, hypoxia would, however, coincide with elevated numbers of quiescent cells, increased cellular debris and local variations of interstitial pressure, all of which are known to impact RT treatment efficacy. It remains to be shown whether resolving only one of these problems, e.g. by increasing the overall level oxygenation, would actually enhance the response of the tumour as a whole. Systems oncology simulations performed for a well-characterized cell line can help to understand the importance of the individual mechanisms through sensitivity analysis better, and to provide a biologically weighted planning tool for the optimization of multimodality treatment scheduling and dosage. 3D tumour spheroids have been shown to recapitulate important features of tumours in vivo, such as gradients in oxygen and nutrients, or 3D cell-cell contact. However, a spheroid model such as the one presented here does not account for the impact of physiology, in particular tumour vasculature, surrounding normal tissue, and the clearance of cellular debris by the immune system. It is thus expected that the translation of this model to the description of tumour response in vivo, will be a complex and challenging undertaking. However, the spheroid model provides an important first step towards this goal conclusions Despite an implementation that deliberately simplified some aspects of biology, the simulation framework presented here was designed to take into account the most important differences between heat-and radiation-induced cell death in 3D tumour spheroids, namely the (in)dependence of cell death on proliferation. The framework was successfully calibrated and validated against experimental data from a commonly used human colon tumour cell line (HCT116). After further development and validation, the framework may provide the basis for more sophisticated simulation of tumour response in vivo, that could be used to optimise fractionation and scheduling of single or multi-modality treatments.

Data availability
The data sets used and/or analysed during the current study are available from the corresponding author on reasonable request. Once the project has been concluded in the group, the entire framework will be uploaded to a public code sharing platform.