Efficient Radial-Shell Model for 3D Tumor Spheroid Dynamics with Radiotherapy

Simple Summary Approximately 50% of patients diagnosed with cancer receive radiotherapy at least once during their disease. Experiments with sophisticated in-cellulo assays to improve radiotherapeutic outcomes are still challenging, and some critical details of tumor cell dynamics still need to be explored. To enhance the informative value of such approaches and support future therapeutic study designs, we developed an efficient mathematical model for three-dimensional multicellular tumor spheroids, which reflect microregions within a large tumor or avascular micrometastases and which are an auspicious experimental framework to pre-assess the curative effect of radio(chemo)therapy. We validate our mathematical model using experimental tumor spheroid growth data of several cell lines with and without radiotherapy and observe equal or better performance than previous models. Moreover, our model allows for efficient parameter calibration within previously reported and/or physiologically reasonable ranges. Based on this data-driven approach, we can explain the mechanism of the characteristic dynamics at small tumor volumes. Abstract Understanding the complex dynamics of tumor growth to develop more efficient therapeutic strategies is one of the most challenging problems in biomedicine. Three-dimensional (3D) tumor spheroids, reflecting avascular microregions within a tumor, are an advanced in vitro model system to assess the curative effect of combinatorial radio(chemo)therapy. Tumor spheroids exhibit particular crucial pathophysiological characteristics such as a radial oxygen gradient that critically affect the sensitivity of the malignant cell population to treatment. However, spheroid experiments remain laborious, and determining long-term radio(chemo)therapy outcomes is challenging. Mathematical models of spheroid dynamics have the potential to enhance the informative value of experimental data, and can support study design; however, they typically face one of two limitations: while non-spatial models are computationally cheap, they lack the spatial resolution to predict oxygen-dependent radioresponse, whereas models that describe spatial cell dynamics are computationally expensive and often heavily parameterized, impeding the required calibration to experimental data. Here, we present an effectively one-dimensional mathematical model based on the cell dynamics within and across radial spheres which fully incorporates the 3D dynamics of tumor spheroids by exploiting their approximate rotational symmetry. We demonstrate that this radial-shell (RS) model reproduces experimental spheroid growth curves of several cell lines with and without radiotherapy, showing equal or better performance than published models such as 3D agent-based models. Notably, the RS model is sufficiently efficient to enable multi-parametric optimization within previously reported and/or physiologically reasonable ranges based on experimental data. Analysis of the model reveals that the characteristic change of dynamics observed in experiments at small spheroid volume originates from the spatial scale of cell interactions. Based on the calibrated parameters, we predict the spheroid volumes at which this behavior should be observable. Finally, we demonstrate how the generic parameterization of the model allows direct parameter transfer to 3D agent-based models.

Simple Summary: Approximately 50% of patients diagnosed with cancer receive radiotherapy at least once during their disease.Experiments with sophisticated in-cellulo assays to improve radiotherapeutic outcomes are still challenging, and some critical details of tumor cell dynamics still need to be explored.To enhance the informative value of such approaches and support future therapeutic study designs, we developed an efficient mathematical model for three-dimensional multicellular tumor spheroids, which reflect microregions within a large tumor or avascular micrometastases and which are an auspicious experimental framework to pre-assess the curative effect of radio(chemo)therapy.We validate our mathematical model using experimental tumor spheroid growth data of several cell lines with and without radiotherapy and observe equal or better performance than previous models.Moreover, our model allows for efficient parameter calibration within previously reported and/or physiologically reasonable ranges.Based on this data-driven approach, we can explain the mechanism of the characteristic dynamics at small tumor volumes.
Abstract: Understanding the complex dynamics of tumor growth to develop more efficient therapeutic strategies is one of the most challenging problems in biomedicine.Three-dimensional (3D) tumor spheroids, reflecting avascular microregions within a tumor, are an advanced in vitro model system to assess the curative effect of combinatorial radio(chemo)therapy.Tumor spheroids exhibit particular crucial pathophysiological characteristics such as a radial oxygen gradient that critically affect the sensitivity of the malignant cell population to treatment.However, spheroid experiments remain laborious, and determining long-term radio(chemo)therapy outcomes is challenging.Mathematical models of spheroid dynamics have the potential to enhance the informative value of experimental data, and can support study design; however, they typically face one of two limitations: while non-spatial models are computationally cheap, they lack the spatial resolution to predict oxygendependent radioresponse, whereas models that describe spatial cell dynamics are computationally expensive and often heavily parameterized, impeding the required calibration to experimental data.Here, we present an effectively one-dimensional mathematical model based on the cell dynamics within and across radial spheres which fully incorporates the 3D dynamics of tumor spheroids by exploiting their approximate rotational symmetry.We demonstrate that this radial-shell (RS) model reproduces experimental spheroid growth curves of several cell lines with and without radiotherapy, showing equal or better performance than published models such as 3D agent-based models.Notably, the RS model is sufficiently efficient to enable multi-parametric optimization within previously reported and/or physiologically reasonable ranges based on experimental data.Analysis of the model

Introduction
The quantification of the biological effects of (radio)therapy is crucial for designing and evaluating treatment protocols and the response prediction for in vivo tumors.Tumor spheroids are the preferential in vitro model to study the possible means of tumor suppression [1][2][3][4][5][6][7], as they reflect many characteristic features affecting tumor growth dynamics, including three-dimensional (3D) reciprocal intercellular and cell-extracellular matrix interactions.Multicellular spheroids are 3D avascular aggregates of several thousand tumor cells that mimic tumor microareas or micrometastases.Due to their 3D structure, spheroids exhibit metabolic gradients of oxygen, nutrients, and waste products, which can strongly modulate the therapy response of cells in addition to 3D cellular interactions [3][4][5][6][7][8].Oxygen deficiency in tumors (hypoxia) is associated with substantial radioresistance [9][10][11][12].The lack of such resistance factors is one reason that two-dimensional in vitro approaches such as the classical clonogenic survival assays reflect the therapeutic response of cancer cells in tissue comparatively poorly [13].Despite the complexity of the incorporated cellular processes, experiments on in vitro spheroids are less laborious and ethically questionable than in vivo animal studies [14][15][16].Spheroid behavior and response to treatment is often evaluated in terms of the development of the spheroid volume or average diameter over time, usually denoted as the growth curve.State-of-the-art long-term spheroid-based assays assess the curative potential of new treatment strategies by classifying each spheroid within a population as either controlled or relapsed due to its individual growth recovery and kinetics, respectively.The spheroid control probabilities and selected spheroid control doses are computed as analytical endpoints by averaging the therapeutic response over ensembles of spheroids for each treatment dose [17][18][19][20], analogous to the tumor control dose employed for in vivo experiments with mice [21,22].Growth curves and spheroid control probabilities depend on the spheroid type (cell line), the size of the spheroids at the onset of treatment, and the applied dose.
While tumor spheroids provide a physiologically more realistic in vitro framework to study tumor growth dynamics and therapeutic response, these experiments are challenging, as typically hundreds of individual spheroids have to be cultured, irradiated with different doses, and subsequently monitored for a sufficient period in order to obtain statistically reliable results for a single treatment arm.Moreover, spheroids are typically composed of several spatial domains resulting from the pathophysiological metabolic gradients, e.g., the oxygen distribution.The most prominent spatial domains developing with increasing spheroid size and volume, respectively, are a secondary necrotic core surrounded by a viable rim.Both are entangled with growth dynamics and therapy response [23,24].However, these inner structures can only be sampled with sufficient accuracy in stained histological sections, which impedes monitoring of the dynamics in individual spheroids over time.
Mathematical modeling can provide access to these dynamics, histomorphological features, and quantitative description of treatment response, and may allow the prediction of the effect of individual therapies [25].Therefore, such models can assist in selecting the most relevant therapeutic parameters, e.g., the range of radiation doses, to design studies and optimize experimental setups.Starting with the phenomenological description of spheroid growth curves, for example by fitting the Gompertz function [26][27][28], several mathematical models have been formulated and simulated with different approaches, complexities, and objectives.Overall, the mathematical models that describe spheroid dynamics can be categorized into nonspatial and spatial models, depending on whether or not they resolve the spatial distribution of cells.Nonspatial models concentrate on free spheroid growth, i.e., spheroid growth without treatment, starting with Greenspan [29], who explored a set of ordinary differential equations (ODE) to describe the dynamics of the radius of the spheroid and its secondary necrotic core over time.Likewise, Grimes et al. [24] used an ODE model to reproduce the free (untreated) spheroid volume growth of different cell lines mainly based on the cellular doubling time and oxygen consumption rate determined in 2D culture.Using a combination of stained histological sections of spheroids and monolayer-based Seahorse extracellular flux analyses, a relationship between the oxygen consumption rate of individual cells and the necrotic radius within the spheroids was derived [23,24].Browning et al. [30] introduced an ODE model to simulate the radius of a spheroid over time.They further added a statistical approach to estimate the effects of expected fluctuations, finding compelling evidence that WM983b spheroids have a limiting size that is independent of the initial seeding density.Alternatively, in order to resolve the spatial structure of spheroids, Brüningk et al. [1] proposed a 3D cellular automaton and reproduced growth experiments without treatment as well as experiments with combined hyperthermia and radiation treatment for a colorectal cancer spheroid model (HCT-116).Amereh et al. [31] used a coupled system of partial differential equations (PDEs) and ODE to simulate the radius as a function of time during the early formation of tumor spheroids without treatment and derived an analytical solution for their model.Ward and King [32] and Jin et al. [33] both used models that explicitly included spatial variations in cell densities and oxygen concentrations based on PDEs to study the phenomenon of growth saturation in tumor spheroids.In both instances, their results were qualitatively compared to experimental data from real-time cell cycle imaging in spheroids [33].Paczkowski et al. [34] employed a nonspatial logistic growth model to describe the growth of homogeneous tumor spheroids and a cellular automaton model representing a two-dimensional (2D) cross-section through a 3D tumor spheroid to demonstrate that prostate cancer spheroids comprising mixed tumor cell populations have different growth kinetics and radiosensitivities than spheroid monocultures.The 2D cross-section was used to reduce the computational cost of a full 3D cellular automaton.
Typically, models of spheroid kinetics face at least one of two challenges.On the one hand, while computationally cheap, nonspatial models lack the resolution of an oxygen gradient across the spheroid required to estimate the oxygen enhancement of radioresistance.On the other hand, while models with genuine 3D resolution can provide metabolic gradients they are computationally expensive, particularly cell-based models that scale with the number of cells.This makes calibration of parameters based on experimental data very challenging, especially as the stochastic nature of cell-based models requires an ensemble of simulation runs for each setup to obtain the average dynamics.
To overcome the challenges of both nonspatial and spatial model types, we propose a radial-shell (RS) model (Figure 1) that is spatially discrete and continuous in time.This model can be seen as a special one-dimensional Markov chain representing the cells' dynamics in and between radial shells.For this, we assume the spheroids to be rotationally symmetric, as suggested by histological sections [18,24] and exploited by several mathematical models before [1,24,[29][30][31]34].We demonstrate that this model is (i) spatially sufficiently resolved to describe the effects of radiotherapy, and (ii) has sufficient computational performance to enable multi-parametric optimization.We show that our model can reproduce free (untreated) growth and radiotherapy response for several cell lines, including estimates of the secondary necrotic core.Its agreement with the experimental data is similar to or better than the models originally proposed for the data.Notably, this is achieved without increasing the number of fit parameters compared to the previous cell-based models.Moreover, almost all model parameters correspond to physiologically meaningful entities and are calibrated within biologically realistic ranges.We further demonstrate that virtually all model parameters can be directly transferred to a corresponding 3D cell-based model, which is potentially valuable for examining stochastic deviations in single-cell behavior as well as rotational symmetry.(a) Sketch of a typical stained histological section of a spheroid with a viable rim (green) surrounding a secondary necrotic core, with cell debris membrane-defect cells (red) suggesting an approximate rotational symmetry of the spheroid.(b) By exploiting the symmetry, the model is able to consider the cell dynamics on radial shells i ∈ [0, 1, 2, ...) and along the distance r from the center of the spheroid at i = 0.The shells have equal radial width dr, and their maximal volume V i is consequently increasing outwards.Each shell i holds a concentration of cells c T (r i ) of type T cells, where T = p denotes proliferation-competent cells (green, proliferate with maximal rate γ, consume oxygen with rate a) and T = n membranedefect cells (red, irreversibly lost competence to proliferate, do not consume oxygen).Note that the concentrations c T ∈ [0, 1] are normalized according to the available space in each shell.(c) Oxygen consumption of proliferation-competent cells leads to a radial decrease of oxygen pressure from the value ρ 0 at the most outer shell at r 0 down to an anoxic threshold ρ an at r an towards the center of the spheroid.Beyond the threshold ρ an , cells are anoxic and become membrane-defect with the anoxic death rate .The volume occupied by membrane-defect cells is reduced with rate δ and all cells are transported inward with rate λ such that the spheroid remains compact.Then, the outer radius R spheroid and necrotic radius R necrotic are defined as the radius of a sphere with a volume equal to the total cell volume and membrane-defect cell volume, respectively.For reference, the shell width dr = κdr * is expressed as multiple κ of a single cell diameter dr * .

Radial-Shell Model (RS Model)
We present a data-driven dynamical model for the growth kinetics of a spheroid that both resolves the radial cell composition of the spheroid and at the same time is efficient enough for manageable parameter optimization; see Figure 1 for an illustration and Appendix C for details.The spheroid is assumed to be rotationally symmetric, as suggested by experimental observations [18,24].The resulting one-dimensional (1D) dynamics are defined by the average concentrations c T (r, t) of type T cells and the partial oxygen pressure, i.e., the oxygen pressure ρ(r, t) along the radius r of the spheroid discretized into radial shells of width dr.
For free (untreated) growth, two cell types T are considered: proliferation-competent cells T = p, and membrane-defect cells T = n.In the model context, proliferationcompetent cells refers to cells with the innate capacity to proliferate with a maximal proliferation rate γ given sufficient oxygen supply and free space in their vicinity.In contrast, membrane-defect cells have undergone secondary necrosis, and as a result have permanently lost their proliferative capacity.These cells do not consume oxygen, and the spheroid volume occupied by their cell bodies and related cell debris declines at a defined rate δ.We further assume as an approximation that proliferation-competent cells consume oxygen according to a constant local oxygen consumption rate a as long as the oxygen pressure is above a (putatively anoxic) threshold ρ an , resulting in an oxygen gradient ρ(r, t) along the radius.At low oxygen pressure ρ(r, t) ≤ ρ an (severe hypoxia), the initially proliferation-competent cells can no longer consume oxygen to maintain their energy homeostasis and proliferation potential.Consequently, they eventually lose their survival capacity and membrane integrity with a defined rate , leading to the emergence of a core of dead cells in the center of the spheroid surrounded by a rim of proliferation-competent cells (see Figure 1 for an illustration).This structure is consistent with the histomorphological observation of a viable rim and a necrotic center in median sections of spheroids with increasing size [3,7,18,24,27].
Notably, not all proliferation-competent cells in the spheroid viable rim actively proliferate simultaneously and in the specific local micromilieu.Subpopulations undergo transient cell-cycle arrest in the 3D cellular environment of the spheroid.However, these cells may re-enter the cell cycle by stimulation via sufficient competence factors, and as such can contribute to relapse.Vice versa, the viable rim may contain membrane-intact and proliferation-incompetent cells that consume oxygen.While these permanently cell-cycle arrested cells, i.e., terminally differentiated or senescent cell populations, can alter the oxygenation within the spheroid and thereby affect the adjacent proliferation-competent cells, they do not directly contribute to the treatment response in terms of control or relapse of a spheroid, as they cannot return to an active cycling phenotype.While it is straightforward to incorporate such cells in our model, we neglect this cell type in order to obtain an effective model with a minimal number of parameters.
A radial shell i ∈ {0, 1, 2, . ..} encompassing all points with radius r ∈ [i 3 ).Note that while the width dr of the shells is constant, their volume increases from the center i = 0 outwards and scales quadratically V i ∼ i 2 dr 3 ∼ r 2 i dr for large radii r i .The average concentration of cell type T in shell i is c T (r i ) ∈ [0, 1] normalized such that a total concentration ∑ T c T (r i ) = 1 represents the maximal dense packing of cells.Thus, c T (r i ) is a relative cell concentration with respect to the available space in the corresponding shell i, and has no unit.The cell concentration added from the proliferation-competent cells T = p in shell i is distributed over shells i − 1, i, and i + 1 proportional to the available free space.Thus, the shell width dr sets the spatial scale for cell transport and the distribution of daughter cells from proliferation.Notably, the width dr is not merely a discretization parameter that can be chosen to be arbitrarily small, as the effective proliferation and consequent growth dynamics depend on dr.Note that proliferation is reduced when the available free space is deficient, analogous to logistic growth.Furthermore, proliferation is assumed to be inhibited when the local oxygen pressure ρ(r, t) is below a hypoxic threshold ρ h ρ an .Further, cells which experience oxygen pressure below the so-called anoxic threshold ρ an become membrane-defect cells with an 'anoxic death rate' .The volume occupied by membrane-defect cells is reduced at a specific rate δ.These transitions reflect diverse biological processes, including various cell death mechanisms ranging from necrosis via necroptosis, ferroptosis, and classical apoptosis, which, for brevity, are not further distinguished in the model.The quantification of these processes requires the critical radius r an at which the oxygen pressure reaches the anoxic threshold ρ(r an ) = ρ an .This radius is determined from the oxygen profile ρ(r), which is obtained as the steady-state solution of the rotationally symmetric diffusion equation for the given distribution of cells c T (r i ), oxygen consumption rate a, and oxygen diffusion constant D ρ .As a boundary condition, we assume a constant oxygen pressure ρ 0 at the surface of the spheroid at r 0 provided by the center of the outermost shell for which the total cell concentration exceeds 0.1; see Appendix D for details.Note that additional metabolic gradients such as nutrition and waste products can be incorporated into the model; for simplicity, we use oxygen here as an effective representative gradient.Finally, to maintain a dense spherical structure, an inward transport with rate λ is implemented to account for the net effect of biomechanical mechanisms present in the spheroid, such as intercellular adhesion and mechanical pressure (see Figure 1 for an illustration).

Dynamics with Radiotherapy
The RS model not only resolves the growth dynamics of an untreated spheroid, it also allows the incorporation of therapeutic effects relevant to overall tumor spheroid response.In particular, in addition to the radial distribution of cells c T (r i ), the model provides the corresponding oxygen profile ρ(r), which is crucial for the radioresistance of cells.Oxygen deficiency in tissues (hypoxia) is a well-known pathological radioresistance factor; oxygen-depleted tumor areas require up to three times higher radiation doses than well-oxygenated tumor areas to achieve the same therapeutic response [10][11][12].Various models exist to describe the effect of radiotherapy on cells [35,36].In general, these models predict the probability S(d eff ) of a cell surviving irradiation with an effective dose d eff , most prominently the linear quadratic (LQ) model S(d eff ) = 1/ exp(α RT d eff + β RT d eff ) [37] with the therapeutic ratio α RT , β RT of the cell line-dependent radiological constants [37][38][39], where the effective dose d eff (d, ρ) ≤ d is considered a function of the actual dose d and local oxygen pressure ρ, reflecting the oxygen enhancement of radioresistance.A prominent choice for the oxygen enhancement ratio (OER) d/d eff is the Alper-Howards-Flanders oxygen equation [9].Mitotic catastrophe has been suggested as a process for radiation-induced cell death in spheroids as a consequence of lethal irreparable DNA damage [40][41][42][43].The incorporation of a wide variety of such mechanisms for the effect of radiation into the RS-model is straightforward.As an example, and to evaluate the performance of the model, we incorporated mitotic catastrophe as recently proposed for a 3D cellular automaton [1].When the spheroid is irradiated with dose d, proliferation-competent cells c p (r i ) are instantaneously converted into damaged cells c d (r i ) with probability 1 − S(d eff (d, ρ(r i ))).Damaged cells proliferate with probability 1 − P mc analogous to proliferation-competent cells, i.e., with rate γ(1 − P mc ); however, they break up with probability P mc during proliferation, i.e., the volume of damaged cells is reduced with a rate γP mc implemented analogously to the reduction of necrotic volume.As previously assumed [1], the probability P mc of mitotic catastrophe takes an initial value P 1 mc < 0.5 and increases to P 2 mc > 0.5 some time after treatment.

Parameter Calibration
The parameters of the RS model can be calibrated based on experimental volume growth curves of spheroids.Further, these parameters can be categorized into four groups: known cell-line independent parameters (D ρ , ρ 0 , ρ an , dr * ), set by fixed values, known cell-line specific parameters (γ, a, α RT , β RT ), calibrated within their reported range (Tables 1 and 2), unknown biological parameters ( , δ, P mc ), and effective parameters (λ, dr) which are fitted within a reasonable range.The reported parameter ranges are taken exclusively from monolayer experiments due to a lack of systematic parameter reporting in spheroid experiments [44]; see Table 1 for a list of references.The specific relationships between the parameters observed in the monolayer experiment and those expected in the spheroid setup are discussed below.
As a simplification for the model, the oxygen diffusion constant D ρ , oxygen pressure ρ 0 at the surface of the spheroid, and anoxic threshold ρ an are set to fixed values of D ρ = 2 × 10 −9 m 2 s −1 , ρ 0 = 100 mmHg, and ρ an = 0 mmHg according to those used in Grimes et al. [23,24] and Brüningk et al. [1].While the oxygen pressure ρ 0 depends in principle on the experimental setup [45,46], we found that within the experimentally observed range of ρ 0 ∈ [80, 120] mmHg the resulting dynamics are insensitive to the exact oxygen pressure ρ 0 (Tables A1 and A2).The proliferation rate γ and the oxygen consumption rate a are fitted within their cell line-specific range as implied from monolayer culture doubling times in the literature, using the latter as an upper bound of the cell cycle time (Table 1).In general, a spheroid's maximal cell proliferation rate γ is considered to be smaller than in the corresponding monolayer culture.However, the observed values are largely dependent on the specifics of the experimental setup [44].Note that for the FaDu cell line we assume a range of ln(2)/γ ∈ [20,40] h, around the reported value of 30 h.In the case of the oxygen consumption rate a, it has been proposed that the value a measured for single cells is consistent with the size of the secondary necrotic core observed in histological sections of spheroids [23,24].As a simple reference for the spatial scale, the average single cell diameter dr * = 16 µm is fixed for all cell lines to an intermediate value of the sizes reported in the literature for the cell lines of interest [17,47,48].The model dynamics are independent of this single cell diameter dr * , as the oxygen consumption rate is defined per cell volume and not per cell number.Neither the anoxic death rate nor the rate δ at which the volume occupied by membrane-defect cells is reduced have been directly measured in experiments for these cell lines, and as such are fitted using a range of several orders of magnitude around the proliferation rate.The inward transport rate λ is fitted in the order of cell transport rates in tissue [49].The shell width dr is simultaneously fitted in the range of several single cell diameters dr = κdr * .Notably, the width dr sets the spatial scale for cell transport and the distribution of daughter cells due to proliferation.In particular, the factor κ can be translated to a spatial neighborhood of a cell; see Section 4.1.
We optimized the mentioned model parameters by minimizing the residuum between the experimentally measured growth curve Rspheroid (t) and the one predicted by the model R spheroid (t) using nonlinear least-squares minimization provided by the Python package lmfit [50].Because the RS model returns the radial cell concentration c T (r i ) for any time point t and for every cell type T, we define the corresponding radius of the spheroid R spheroid via the total cell volume spheroid .To take into account the inner dynamics of the spheroid, we consider the outer radius Rspheroid and the secondary necrotic core simultaneously, minimizing the residuum between the radii of the secondary necrotic core based on spheroid growth experiments Rnecrotic and the model R necrotic .While R necrotic is defined analogous to R spheroid for the model, based on the volume of membranedefect cells T = n, measuring the necrotic radius Rspheroid requires histological sections, which impedes monitoring those dynamics over time in individual spheroids.A method for estimating the oxygen consumption rate a from the outer radius Rspheroid and necrotic radius Rnecrotic has been devised previously (Equation (A42)), and these estimates seem to be consistent with extracellular flux measurements of the oxygen consumption rate [23,24]; thus, we reverse this method to infer the corresponding necrotic radii Rnecrotic (t) from the growth curve Rspheroid (t) using published values for the consumption rate a. Notably, while the radii of the spheroid R spheroid and R necrotic are similar in general, they are not equal to the radii r 0 and r an from the oxygen profile.We define t = 0 d as the first day of spheroid monitoring, when the spheroid has already formed.However, it needs to be emphasized that experiments usually start with the seeding of cells a few days prior.The seeding procedure and forming of the spheroid may differ between experiments; these details are neglected in our model [31].Finally, the initial condition for the model reflecting an in vitro spheroid with Rspheroid (t = 0) is set by the following procedure: the radial distribution of cell concentrations c p (r i ) is initialized with a generic sigmoid function corresponding to a spheroid with much smaller radius, e.g., reflecting only 20% of the initial experimental spheroid volume, such that there is no anoxic core; subsequently, the model dynamics are integrated until the simulated spheroid reaches the size of the experimental one at t = 0, and the corresponding distribution of cell concentrations is used as initial condition (see Appendix J for details).
Subsequently, the parameters α RT , β RT , and P mc related to radiotherapy are calibrated: the cell-specific radiosensitivities α RT and β RT are taken from the ranges previously reported from clonogenic survival assays [37][38][39].The initial and final probabilities of mitotic catastrophe P 1 mc and P 2 mc are optimized based on the dynamics of the outer radius observed at high doses, e.g., 10 Gy for HCT-116 and 20 Gy for FaDu, for which spheroids do not relapse anymore.This approach is analogous to the calibration of the other model parameters, except that only P 1 mc and P 2 mc are optimized while the other parameters are fixed to the values calibrated based on untreated growth.
For the purposes of comparison, the parameters used in the models proposed along with the experimental data are listed in Tables 1 and 2. For the cell-based model of HCT-116, only the probabilities of necrosis • dt and the reduction of the volume occupied by membrane-defect cells δ • dt are reported [1].For a Monte Carlo step size of dt = 1 h, this translates to the corresponding rates displayed in Table 1.In the nonspatial model [24], the whole secondary necrotic core is removed within a time step dt = ln(2)/γ, from which we estimate a lower bound for the product • δ ≥ γ.
Table 1.The RS model reproduces experimental growth curves with parameters in a realistic biological range.For each experiment, the estimated parameter values and fit result of the calibrated RS Model, literature values from monolayer experiments for the parameters, and fit results from previous models (except for FaDu, for which no model has been proposed) are displayed.The parameters (left to right) are: doubling time ln(2)/γ with the proliferation rate γ, oxygen consumption rate a, anoxic death rate , rate δ at which the volume occupied by membrane-defect cells is reduced, the width of the radial shells dr = κdr * expressed in multiples κ of the cell diameter dr * , and the inward transport rate λ, displayed as the velocity λ • dr.The goodness of fit is assessed by the coefficient of determination R 2 for the volume V spheroid .Note that the two sets of lines on HCT-116 refer to separate experiments which exhibit different spheroid growth curves, probably due to discrepancies in the experimental setup.For the FaDu cell line, we assume a range of ln(2)/γ ∈ [20,40] h, which is around the reported value of 30 h. ∼ 10.6

Cell Line Parameter Estimation ln(2)/γ [h] a [mmHg/s] [1/h]
Table 2.The RS model reproduces the experimentally observed spheroid dynamics after radiation for HCT-116 and FaDu.For HCT-116, the same radiosensitivities as in Brüningk et al. [1] were chosen and the probabilities of mitotic catastrophe P mc were optimized based on the highest dose of 10 Gy while adopting the other parameters for untreated growth (Table 1).Similarly, for FaDu the probabilities of mitotic catastrophe were optimized for 20 Gy, while apt radiosensitivities were chosen from the range observed in previous clonogenic survival assays [17,59] and optimized for 5 Gy.The goodness of fit was assessed using the coefficient of determination R 2 for the volume V spheroid .The bold R 2 values indicate the experiments used for validation (not included in the calibration).[1] (black crosses) even better than the model proposed along with the data (light green dashed).Note that the necrotic radius during the experiment can be estimated (gray crosses; see Section 3 for details) and that the RS model additionally reproduces this necrotic radius (red).(c) Volumes corresponding to the radii displayed in (a) in a semi-log plot.Multiples of 2 and 5 of the starting volume V t=0 are highlighted with horizontal gray dashed and gray dotted lines, respectively.The RS model predicts the cell concentration over the radius at any time point.Examples show the cell distributions for (b) the first time point t = 0 d and (d) the last day t = 21 d of the experiment.Displayed are the cell concentrations of proliferationcompetent cells c p (green solid), membrane-defect cells c n (red solid), total concentration of cells (black dotted), and oxygen profile (blue solid, right y axis).The marker at the bottom of the panels (b,c) indicates the necrotic and outer radius estimated from the model (red and green crosses) and from the actual experimental data (gray and black arrows), respectively.The used parameters are reported in Table 1.

Results
In order to demonstrate the feasibility and performance of the RS model, we first calibrated the model for five different cancer cell lines (colorectal HCT-116 and LS-174T, breast MDA-MB-468, squamous cell SCC-25, and pharynx carcinoma FaDu) from six individual experiments of untreated growth [1,18,24] and compared the goodness of fit with the models previously proposed for these experiments (Table 1).The experiments included the cell lines HCT-116, LS 174T, MDA-MB-468, and SCC-25, which have been previously compared to a nonspatial model [24], HCT-116 with a comparison to a cellular automaton model [1], and FaDu, for which no corresponding model has been proposed [18].The goodness of fit was quantified by the coefficient of determination R 2 with respect to the spheroid volume, as reported for previous models [24].In all cases of untreated growth, the RS model yields similar or higher goodness of fit than the corresponding models proposed for these experimental data.As an example, experimental spheroid radii Rspheroid (t) and model predictions R spheroid (t) are displayed for HCT-116 in Figure 2a, while the corresponding figures for the other cell lines can be found in Figures A1-A5 along with the same results reported in terms of spheroid volume.Furthermore, while the goodness of fit is only reported with regard to the spheroid radius, the RS model is additionally calibrated to reproduce the dynamics of the secondary necrotic core (Figure 2a).Note that in the RS model the secondary necrotic core is a feature emerging from the dynamics within the spheroid, while for the nonspatial model proposed for the other cell lines [24] the necrotic radius R necrotic is simply identified with the anoxic radius r an .In particular, the distribution of membrane-defect and proliferation-competent cells over the radius of the spheroid is resolved, as shown by the examples of the first and last time point in Figure 2b,d; see Video S1 for the full dynamics over time.
Subsequently, the model was applied to experimental results after radiotherapy, available for HCT-116 at doses of 2, 5, 10 Gy [1] and for FaDu at doses of 2.5, 5, 20 Gy [18]; see Table 2 for the parameters and goodness of fit analogous to Table 1.For HCT-116, the RS model reproduces the dynamics of the spheroid after treatment as well as the model proposed along with the data.Note that only the probabilities for mitotic catastrophe P mc are additionally calibrated based on the growth curve at 10 Gy (Figure 3), while the data for 2 and 5 Gy validate the model prediction (Figures 4 and 5).The RS model reproduces the experimental results for FaDu with goodness of fits comparable to the case of HCT-116 (Figures A6-A8).Because a very wide range of radiosensitivities α RT , β RT has been reported in the literature for FaDu [17,59], we tried to choose apt values with respect to the dynamics at 5 Gy, which is the dose most sensitive to these parameters.
We highlight that the spatial resolution and agreement with the experimental results is achieved without increasing the number of fit parameters compared to the cell-based model.All the model parameters are calibrated within previously reported or physiologically reasonable ranges; see the comparison with the values reported in the literature in Table 1.Considering the challenge of parameter identifiability typical of data-driven models, the calibrated parameters reveal additional insights into cell-specific characteristic of spheroid dynamics.From the analytic solution of an ODE model of the effective spheroid dynamic, it can be inferred that the slope of a typical growth curve is determined by the product of the proliferation rate γ and shell width dr = κdr * , which sets the spatial scale for cell transport and distribution of daughter cells from proliferation (see Appendix C and Figure A11).This means that at intermediate-to-large spheroid sizes R spheroid > 3 • dr, for which effectively only an outer shell of the spheroid is able to proliferate, the proliferation rate (γ in the case of untreated growth and γ(1 − 2P mc ) after radiation) and dr cannot be inferred individually from the growth curves, only from their mutual product.However, the radial dynamics are exponential at a small spheroid radius R spheroid ∼ dr.The exponential phase corresponds to the scenario in which all proliferation-competent cells within the spheroid can proliferate and the exponent is essentially set by the proliferation rate.This behavior is particularly pronounced for FaDu after strong radiation for R spheroid < 100 µm (Figure A6), and to some degree for MDA-MB-468 (Figure A2.This not only allows us to determine dr = κdr * and γ independently for these cases, it demonstrates the biological relevance of the spatial scale dr.In particular, for the two cases of MDA-MB-468 and LS-174T with relatively large κ, there is a lower bound for κ ≤ 5.4 and ≤ 5.7, respectively, assuming the largest proliferation rate within the range observed in monolayer, i.e., ln(2)/γ > 30 h for both cell lines.Thus, we predict that for these two cell lines the exponential growth regime should be observable for spheroid sizes in the range of R spheroid ∈ [dr, 3dr] ≈ [90, 260] µm, which is currently slightly below the experimentally observed sizes.Figure A12 [1] (black crosses) as effectively as the model proposed along with the data (light green dashed).We used the parameters calibrated for untreated growth, only fitted the probabilities of mitotic catastrophe P 1 mc , P 2 mc , and adopted the radiosensitivities α RT and β RT from Brüningk et al. [1].Note that for this regime no secondary necrotic core is expected (experimental estimate shown as gray crosses), which is reproduced by the RS model (red).The parameters we used can be found in Tables 1 and 2 [1] (black crosses) even better than the model proposed along with the data (light green dashed).We used the parameters calibrated for untreated growth, the values for P 1 mc and P 2 mc fitted at 10 Gy, and the values α RT and β RT from Brüningk et al. [1].Note that the necrotic radius during the experiment can be estimated (gray crosses) (see Section 3 for details) and that the RS model reproduces this necrotic radius (red).Additionally, note that the necrotic radius of the RS model appears jagged due to the spatial discretization of r 0 .The used parameters are reported in Tables 1 and 2 .We used the parameters calibrated for untreated growth, the values for P 1 mc and P 2 mc , and the model proposed along with the data (light green dashed) at 10 Gy, as well as the values α RT and β RT from Brüningk et al. [1].The necrotic radius during the experiment can be estimated (gray crosses) (see Section 3 for details) and the RS model reproduces this necrotic radius (red).Note that the necrotic radius of the RS model appears jagged due to the spatial discretization of r 0 .The used parameters can be found in Tables 1 and 2. (c) Volumes corresponding to the radii displayed in (a) in a semi-log plot.Multiples of 2 and 5 of the starting volume V t=0 are highlighted with horizontal gray dashed and gray dotted lines, respectively.Additionally, the RS model predicts the cell concentration over the radius at any time point.Examples show the cell distribution for (b) the first time point t = 0 d and (d) the last day t = 21 d of the experiment.Displayed are the cell concentrations of proliferation-competent cells c p (green solid) and membrane-defect cells c n (red solid), the total concentration of cells (black dotted), and the oxygen profile (blue solid, right y axis).The marker at the bottom of the panels (b,c) indicates the necrotic and outer radius estimated from the model (red and green crosses) and from the actual experimental data (gray and black arrows), respectively.The used parameters are reported in Table 1.

Parameter Transfer to a Cellular Automaton
The parameters of the RS model can be translated to a cell-based model which reflects analogous cell mechanisms, as the processes incorporated in the RS model are generic.i.e., cells consume oxygen with a rate a and proliferate with a rate γ given sufficient oxygen and free space in their neighborhood.The size of the directionally unbiased interaction neighborhood considered for placing potential offspring of a cell in the cellular automaton is directly related to the shell width of the RS model.In addition, cells with insufficient oxygen supply become membrane-defect cells at a rate , the volume occupied by these membranedefect cells is reduced with a rate δ, and all cells move inwards to keep the spheroid composition dense.A cell-based model calibrated in this way displays similar dynamics on average as the RS model while additionally exhibiting statistical variation originating from the stochastic nature of the cell-based processes.This can be used to theoretically predict spheroid control probabilities and spheroid control doses for irradiation therapies, where spheroid relapse and control are distinguished and can serve as as basis to optimize experimental design for in vitro dose ranges and potentially fractionated radiation.
To demonstrate parameter transfer, we chose a cellular automaton as an example, similar to that of Brüningk et al. [1] (see Appendix H for details) and used the FaDu cell line.For brevity, we focused on untreated growth.As explained above, almost all parameters of the RS model can be directly transferred to the cellular automaton; the only exceptions are the inwards transport rate λ and shell width dr = κdr * , which are set through multiples κ of a single cell diameter dr * .For technical reasons involving computational speed (see Appendix G for details), we chose to replace the inwards transport rate with a rate-independent inwards shuffle leading to the requirement that the Monte Carlo time step length dt of the cellular automaton be sufficiently small; see Appendix H for details.
Additionally, the shell width dr in the RS model, and consequently κ, defines the size of the interaction neighboorhood with the mean distance d neigh in the 3D cellular automaton, i.e., it determines the maximal distance at which a cell's daughter cell can be placed.In contrast to classical cellular automata, we allow higher orders of the typical Moore and von Neumann neighboorhoods and random alternations between them, where the order k is a suitably chosen integer.We demonstrate that there is an almost perfect linear dependency between κ and the mean cell distance d neigh of a chosen interaction neighborhood κ ≈ 1.19d neigh − 0.29; see Figure A14 and Table A3 for more details.We use this as a basis to specify an interaction neighborhood for the example of the FaDucell line; a neighborhood which alternates 50:50 between a Moore neighboorhood of order k = 3 (d neigh = 3.34) and k = 4 (d neigh = 4.31) corresponds to κ ≈ 4.26, which is close to the value of κ = 4.31 calibrated with the RS model for FaDu.We find good agreement for the radii R spheroid and R necrotic over time between both models as well as for the cell concentrations (Figure 6).

Discussion
We propose an effectively 1D radial-shell (RS) model for simulating 3D tumor spheroid dynamics with and without radiotherapy while taking into account the oxygen-dependent radiosensitivity of cells.By deliberately exploiting the approximate radial symmetry of spheroids as evident in the data, we are able to considerably reduce the computational complexity of our model without compromising the quality of agreement between model and data compared to previous genuinely 3D cell-based models.We show that the proposed RS model reproduces untreated growth and radiotherapy response experiments for many different cell types and doses, including estimates of the secondary necrotic core, while restricting the parameters to biologically realistic ranges.At the same time, our model has sufficient computational performance to allow for multi-parametric optimization, enabling fits which are similar or better than those obtained from models originally proposed for the data.Additionally, we have shown through examples that there is a parameter mapping between the new RS model and related 3D cell-based models that leads to full agreement of the two models with respect to common experimentally relevant observables.In this way, we open up a way for efficient multi-parametric optimization of those cell-based models while at the same time demonstrating the adequacy of dimensional reduction.
The RS model is an effective implementation of key cell-based mechanisms underlying tumor spheroid dynamics and oxygen-dependent therapy response that is parameterized according to biological interpretables and to some extent experimentally identifiable quantities.In particular, we show that the shell size dr = κdr * in the RS model translates to the size of the neighborhood in a 3D cellular automaton, and as such can be biologically interpreted as an effective range of proliferation.To guide experimental identification, we analytically predict the influence of the effective range dr on different spheroid growth regimes in the model; exponential growth for small radii is superseded by linear growth for larger radii, and the position of the area of transition is essentially determined by the effective range of proliferation (or equivalently, the shell size).Further, the increase in the exponential phase is set by the proliferation rate alone, while the gradient of the linear growth is determined by the product of proliferation rate and effective range of proliferation, which means that parameter identification is enhanced if observation of at least a few small radii can be included in the calibration process.The effect of the range of proliferation dr is observable in the experimental data for FaDu spheroids after strong irradiation.While the spheroids in the available experiments for the other considered cell lines are too large to exhibit the effect of dr, we predict that this effect should be experimentally observable for MDA-MB-468 and LS-174T.Based on our calibrations for these two cell lines, we estimate a range of spheroid diameters 180-520 µm in which the transition between linear and exponential growth regime is expected, which could be validated by future experiments.
It is known that there are several other cell states and intercellular processes present in a spheroid beyond those incorporated into the presented RS model, including proliferationincompetent yet oxygen consuming membrane-intact cells, additional metabolic gradients including waste products, cell cycle arrest, and complex pathways for cell death and destruction due to anoxia or radiation.However, our goal is not to reflect all known biological processes within a spheroid but to propose an effective model which, as demonstrated, is sufficient to reproduce the available experimental growth curves based on a few key mechanisms.The model can be easily extended with additional cellular mechanisms if systematic data related to these mechanisms are available.Moreover, it can be used as a framework to include and theoretically analyze the potential effects of additional determinants, such as genetic heterogeneity of the cells modifying their motility, growth, and radioresponse characteristics, cellular plasticity [60,61], stochasticity in cellular response behavior, and loco-regional and spatio-temporal pathophysiological phenomena of tumor tissues such as chronic and acute hypoxia.While it is straightforward to incorporate additional biological processes into the RS model (e.g., the constant outer oxygen pressure ρ 0 may be replaced by a time-dependent function ρ 0 (t) to mimic acute and intermittent hypoxia), these extensions require the introduction of additional parameters.This is only sensible if these parame-ters are already known from independent experiments or if their impact is observable in complementary measurements which can be compared to simulation results, e.g., stained histological sections [23,24,46,62], cell counts [46,62,63], cell cycle imaging [33], or genetic heterogeneity.Furthermore, while estimates of several cell line-specific parameters can be found in the literature, e.g., for proliferation rates ( [51][52][53][54][55][56][57] and DSMZ (www.dsmz.de)),most of these originate from monolayer experiments.Moreover, the effective values of such parameters in spheroid experiments depend not only on the cell line but on the details of the experimental setup, e.g., oxygenation, nutrition concentration, and seeding procedure.This becomes evident when comparing the two growth experiments without treatment on HCT-116 [1,24], which display considerably different spheroid growth curves and consequently result in different estimated proliferation rates for the RS model.This highlights the importance of systematic high quality spheroid experiments with detailed parameter reporting [44].
In the future, we plan to apply and verify the RS model in several scenarios.First, the parameters can be calibrated based on spheroids of a particular initial size, with the calibrated model consequently used to predict the untreated growth and therapy response of larger or smaller spheroids.This can extend the informative value of spheroid experiments and represent a first step in statistical projection from spheroid control probabilities to tumor control probabilities, as spheroids may be considered avascular microareas of varying size within a macroscopic tumor.Second, the parameters resulting from calibration of the RS model can help to uncover potential differences between spheroid types and experimental conditions as a way to stimulate extended studies and guide future experimental design.For instance, the anoxic death rate, the rate at which the volume occupied by membrane-defect cells is reduced, and the inwards velocity are all parameters of the RS model that can be readily interpreted in biological terms but which have not yet been directly identified experimentally.Third, the parameters can be calibrated based on spheroids treated with a single dose, and the calibrated model can then be used to predict the response to fractionated treatment.While fractionated treatment is more realistic than single-dose treatment, the corresponding in vitro experiments with spheroids are more laborious.In this context, model simulations can help to identify the relevant dose range, thereby reducing the required number of experiments.Finally, our results suggest that the proposed model has the potential to guide the design of future experiments at experimentally unseen doses.We have demonstrated for two examples that, following calibration based on curves from untreated growth and at high radiation dose, the model is able to predict the therapy response at intermediate doses for given radiosensitivities α RT and β RT from monolayer experiments.A prerequisite for this application is the determination of effective radiosensitivities in spheroids and their relation to the corresponding monolayer values.While it should be possible to fit such effective radiosensitivities based on at least two growth curves at intermediate doses, a more promising route would be to independently extract these values from spheroid control probabilities using a suitable statistical model.
The proposed RS model framework stands out through its computational efficacy while maintaining high conformance to genuine 3D cell-based models.It can serve as an expedient means to theoretically analyze and question experimental data by providing estimates of important influencing factors and allowing extrapolation to experimentally unexamined conditions.In particular, it could be used in the future to optimize fractionation and scheduling of single or multi-modality treatments as well as to explore the consequences of competing hypotheses on the cellular therapy response in order to gain understanding and guide future experiments.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/cancers15235645/s1,Supporting Information Movie S1 provides an example illustration of the RS model for untreated growth experiment HCT-116 (MP4).The RS model reproduces the outer radius (dark green) for untreated growth of the cell line HCT-116, with experimental data from Grimes et al. [24] (black crosses) even better than the model proposed along with the data (light green dashed).Note that the necrotic radius during the experiment can be estimated (gray crosses; see main text Section 3 for details) and that the RS model additionally reproduces this necrotic radius (red).(c) Volumes corresponding to the radii displayed in (a) in a semi-log plot.Multiples of 2 and 5 of the starting volume V t=0 are highlighted with horizontal gray dashed and gray dotted lines, respectively.Additionally, the RS model predicts the cell concentration over the radius at any time point.Examples show the cell distribution for (b) the first time point t = 0 d and (d) the last day t = 17 d of the experiment.Displayed are the cell concentrations of proliferation-competent cells c p (green solid) and membrane-defect cells c n (red solid), the total concentration of cells (black dotted), and the oxygen profile (blue solid, right y axis).The marker at the bottom of the panels (b,c) indicates the necrotic and outer radius estimated from the model (red and green crosses) and from the actual experimental data (gray and black arrows), respectively.The used parameters are reported in Table 1.[24] (black crosses) even better than the model proposed along with the data (light green dashed).Note that the necrotic radius during the experiment can be estimated (gray crosses; see main text Section 3 for details) and that the RS model additionally reproduces this necrotic radius (red).(c) Volumes corresponding to the radii displayed in (a) in a semi-log plot.Multiples of 2 and 5 of the starting volume V t=0 are highlighted with horizontal gray dashed and gray dotted lines, respectively.Additionally, the RS model predicts the cell concentration over the radius at any time point.Examples show the cell distribution for (b) the first time point t = 0 d and (d) the last day t = 12 d of the experiment.Displayed are the cell concentrations of proliferation-competent cells c p (green solid) and membrane-defect cells c n (red solid), the total concentration of cells (black dotted), and the oxygen profile (blue solid, right y axis).The marker at the bottom of the panels (b,c) indicates the necrotic and outer radius estimated from the model (red and green crosses) and from the actual experimental data (gray and black arrows), respectively.The used parameters are reported in Table 1.[24] (black crosses) even better than the model proposed along with the data (light green dashed).Note that the necrotic radius during the experiment can be estimated (gray crosses; see main text Section 3 for details) and that the RS model additionally reproduces this necrotic radius (red).(c) Volumes corresponding to the radii displayed in (a) in a semi-log plot.Multiples of 2 and 5 of the starting volume V t=0 are highlighted with horizontal gray dashed and gray dotted lines, respectively.Additionally, the RS model predicts the cell concentration over the radius at any time point.Examples show the cell distribution for (b) the first time point t = 0 d and (d) the last day t = 7 d of the experiment.Displayed are the cell concentrations of proliferation-competent cells c p (green solid) ans membrane-defect cells c n (red solid), the total concentration of cells (black dotted), and the oxygen profile (blue solid, right y axis).The marker at the bottom of the panels (b,c) indicates the necrotic and outer radius estimated from the model (red and green crosses) and from the actual experimental data (gray and black arrows), respectively.The used parameters are reported in Table 1.The marker at the bottom of the panels (b,c) indicates the necrotic and outer radius estimated from the model (red and green crosses) and from the actual experimental data (gray and black arrows), respectively.The used parameters are reported in Table 1.The marker at the bottom of the panels (b,c) indicates the necrotic and outer radius estimated from the model (red and green crosses) and from the actual experimental data (gray and black arrows), respectively.The used parameters are reported in Table 1.[18] (black crosses) as effectively as the model proposed along with the data (light green dashed).We used the parameters calibrated for untreated growth, only fitted the probabilities of mitotic catastrophe P 1 mc , P 2 mc , and adopted the radiosensitivities α RT and β RT from the ranges observed in Chen et al. [17] and Xu et al. [59].Note that no necrotic core is expected for this regime (experimental estimate shown as gray crosses), which is reproduced by the RS model (red).The used parameters can be found in Tables 1 and 2. (c) Volumes corresponding to the radii are displayed in (a) in a semi-log plot.Multiples of 2 and 5 of the starting volume V t=0 are highlighted with horizontal gray dashed and gray dotted lines, respectively.Additionally, the RS model predicts the cell concentration over the radius at any time point.Examples show the cell distribution for (b) the first time point t = 0 d and (d) the last day t = 33 d of the experiment.Displayed are the cell concentrations of proliferation competent cells c p (green solid) and membranedefect cells c n (red solid), the total concentration of cells (black dotted), and the oxygen profile (blue solid, right y axis).The marker at the bottom of the panels (b,c) indicates the necrotic and outer radius estimated from the model (red and green crosses) and from the actual experimental data (gray and black arrows), respectively.The RS model reproduces the outer radius (dark green) for the FaDu cell line with 5 Gy radiation, experimental data from Chen et al. [18] (black crosses).We used the parameters calibrated for untreated growth, the values for P 1 mc and P 2 mc fitted at 20 Gy, and adopted the radiosensitivities α RT and β RT from the ranges observed in Chen et al. [17] and Xu et al. [59].Note that no necrotic core is expected for this regime (experimental estimate shown as gray crosses), which is reproduced by the RS model (red).The used parameters can be found in Tables 1 and 2. (c) Volumes corresponding to the radii are displayed in (a) in a semi-log plot.Multiples of 2 and 5 of the starting volume V t=0 are highlighted with horizontal gray dashed and gray dotted lines, respectively.Additionally, the RS model predicts the cell concentration over the radius at any time point.Examples show the cell distribution for (b) the first time point t = 0 d and (d) the last day t = 21 d of the experiment.Displayed are the cell concentrations of proliferation competent cells c p (green solid) and membranedefect cells c n (red solid), the total concentration of cells (black dotted), and the oxygen profile (blue solid, right y axis).The marker at the bottom of the panels (b,c) indicates the necrotic and outer radius estimated from the model (red and green crosses) and from the actual experimental data (gray and black arrows), respectively.The RS model reproduces the outer radius (dark green) for the FaDu cell line, with 2.5 Gy radiation experimental data from Chen et al. [18] (black crosses).We used the parameters calibrated for untreated growth, the values P 1 mc and P 2 mc fitted at 20 Gy, and adopted the radiosensitivities α RT and β RT from the ranges observed in Chen et al. [17] and Xu et al. [59].The necrotic radius of the experiments can be estimated (gray crosses; see Section 3 for details).The RS model can predict the necrotic radius (red) in a good manner.Note that the necrotic radius of the RS model appears jagged due to the spatial discretization of r 0 .The used parameters can be found in Tables 1 and 2. (c) Shows the same data plotted for the volume in a semi-log plot.Multiples of 2 and 5 of the starting volume V t=0 are highlighted with horizontal gray dashed and gray dotted lines, respectively.Additionally, the RS model predicts the cell concentration over the radius at any time point.Examples show the cell distribution for (b) the first time point t = 0 d and (d) the last day t = 14 d of the experiment.Displayed are the cell concentrations of proliferation competent cells c p (green solid) and membrane-defect cells c n (red solid), the total concentration of cells (black dotted), and the oxygen profile (blue solid, right y axis).The marker at the bottom of the panels (b,c) indicates the necrotic and outer radius estimated from the model (red and green crosses) and from the actual experimental data (gray and black arrows), respectively.

Appendix C. Equations for the RS Model
The 1D radial-shell (RS) model describes the dynamics of concentrations c T (r i ) of cells of different types T across the radial shells of a spheroid, where the i-th shell is centered at r i = (i + 1/2)dr and has width dr.In accordance with the dynamics of individual cells, the dynamics of the concentrations are governed by the following processes, some of which depend on the local oxygen concentration ρ(r i ): • Net gain by proliferation of proliferation-competent cells, T = p, and of damaged cells, T = d, with probability 1 − P mc : f γ (r i , T, ρ) • Net gain due to inward transport for all cell types: f λ (r i , T) • Anoxic death of proliferation-competent and damaged cells T ∈ {p, d}: f (r i , T, ρ) • Reduction of volume occupied by membrane-defect cells: f δ (r i ) • Death due to mitotic catastrophe for damaged cells, T = d, with probability P mc : For the above, the time scales are set by corresponding rates, that is, the proliferation rate γ, inward transport rate λ, anoxic death rate , and reduction rate δ.The time-continuous dynamics of the RS model are then provided by a set of ordinary differential equations: with the individual terms detailed as follows.First, the net gain due to proliferation f γ describes the increase γc T (r i )V i of proliferation-competent cells T = p due to proliferation with rate γ proportional to the number of cells in the shell c T (r i )V i .This proliferation is limited by the available free space in the immediate neighborhood of the cells, which is incorporated using the term l ogistic (i): , where (A3) This means that the proliferation from shell i is linearly reduced following l ogistic (i) < 1 whenever the total available free space in shell i and in its directly neighboring shells i + 1 and i − 1 is less than the volume V i of the shell.Furthermore, the added cell volume from proliferation in shell i is distributed across shells i − 1, i, i + 1 proportionally to the available space in the respective shells.The fraction of added cell volume originating from proliferation in shell i o that a neighboring target shell i τ ∈ {i o − 1, i o , i o + 1} receives is expressed by the term n eighbouring (i o , i τ ): This fraction is determined by the ratio of the free volume of the target shell i τ and the total free volume in all potential origin shells {i o − 1, i o , i o + 1}.Consequently, the net gain for a shell i o encompasses not only the gain from its own proliferation but from the proliferation of the neighboring shells, as incorporated by the sum ∑ 1 k=−1 c T (r i+k ) in Equation (A2).This gain of cell volume in shell i is translated to a corresponding cell concentration in by the normalization 1/V i .The special cases i o ∈ {−1, 0} in Equation (A5) numerically ensure a von Neumann boundary condition at the center of the spheroid.Finally, proliferation is assumed to only be possible if the local oxygen pressure is above the hypoxic threshold ρ h , which is implemented in Equation (A2) by the term 1 − L r h (r i ) with (A6) This term restricts proliferation to cells above the hypoxic radius r h at which ρ(r h ) = ρ h .Note that this does not change the cell type of proliferation-competent cell T = p; it only switches off their proliferation when they experience hypoxic conditions ρ an < ρ i < ρ h .If the the hypoxic radius r h falls into shell i, r h ∈ [r i − dr 2 , r i + dr 2 ], proliferation in this shell is proportional to the fraction of the non-hypoxic volume 1 − All shells located more centrally have no proliferation due to hypoxia, while the proliferation in all shells located more outwardly is unaffected by hypoxia.Note that we adopted the small value ρ h = 0.5 mmHg from Grimes et al. [23,24] and obtained virtually the same results as for ρ h = ρ an = 0 mmHg, which we consequently used in the main text.
Second, the term P mc f mc in Equation (A1) reflects cell death due to mitotic catastrophe on the part of cells T = d damaged by irradiation with This means, that damaged cells T = d are destroyed and immediately removed from the corresponding shell with probability P mc during proliferation.The term f mc reflects the cells which try to proliferate; refer to the net gain due to proliferation in Equation (A2).In turn, the term (1 − P mc ) f γ (r i , d, ρ) represents the part of the damaged cells which proliferate normally.
Note that there is effective outward transport due to proliferation, where the added cell volume is distributed over all neighboring shells (Equation (A2)).Because the shell volume increases with increasing radius and the distribution is based on the available free volume (Equation (A5)), this inward transport is on average directed outwards.Thus, an effective inward transport f λ (r i , T) is incorporated in Equation (A1) including an inward transport rate λ, which reflects all biological processes, e.g., cell adhesion, that keep the tumor compact despite the reduction of volume occupied by membrane-defect cells in the center and proliferation at the surface.The net gain from the inward transport f λ contains two parts, namely, the cell volume migrating from shell i + 1 to shell i and the cell volume migrating from shell i to shell i − 1.For the first part, the cell concentration of the more outer shell i + 1 is translated to a concentration in shell i via the ratio of the volumes V i and V i+1 .Further, the inward transport is linearly restricted by depending on the free volume in shell i, 1 − ∑ T c T (r i ).
Finally, the term f (r i , T, ρ) expresses the death of proliferation-competent and damaged cells T ∈ {p, d} due to anoxia: converting them into membrane-defect cells T = n at a rate .This process is restricted to regions with insufficient oxygen, i.e., to radii below the anoxic radius r an with ρ(r an ) = ρ an , which is incorporated by the factor L r an (r i ) (Equation (A6)) in a manner analogous to the effect of hypoxia during proliferation.In turn, the volume occupied by membrane-defect cells T = n is reduced with rate δ according to f δ in Equation (A1) In order to have no flow of cell volume at the center r = 0 and no inflow of cell volume at the outermost radius r = r max , boundary conditions have to be enforced for the dynamics.Numerically, these boundary conditions are implemented via ghost shells beyond these boundaries at r −2 , r −1 , r n+1 , r n+2 , which have constant cell concentrations and volume.At the center of the tumor, i ∈ {−2, −1}, we set ∑ T c T (r −1 ) = ∑ T c T (r −2 ) = 1 and V −1 = V −2 = V i=0 and use special cases in n eighbouring and l ogistic (Equations (A3) and (A5)).At the outer edge r max , we set a large volume V n+1 = V n+2 = 10 8 and zero cell concentration ∑ T c T (r n+1 ) = ∑ T c T (r n+2 ) = 0 such that cell volume is only transported outside.Note that r max has to be chosen as sufficiently large to ensure ∑ T c T (r max ) ≈ 0 and avoid an impact of the outer boundary on the simulation.
The Python package solve_ivp from scipy.integrate was used to integrate the RS model equations.We set the integration methods to an Explicit Runge-Kutta method of order 3.
for r > r an , with the boundary conditions ρ(r 0 ) = ρ 0 , ∂ρ(r an ) ∂r = 0, ρ(r an ) = ρ an , (A17) These boundary condition specify that at the tumor surface at radius r 0 the oxygen pressure equals ρ 0 .Additionally, the oxygen profile converges at radius r an to the anoxic threshold ρ an , below which cells cannot consume oxygen any more.We integrated Equation (A16) using a computer algebra program: Due to the complexity of Equation (A16), we only used two boundary conditions: ρ(r 0 ) = ρ 0 and ρ(r an ) = ρ an .To further solve the resulting Equation (A18), we interpret the cell concentrations c T (r i ), which represent the average cell concentration in the discrete radial shells with width dr centered at r i = (0.5 + i)dr, as space-continuous step functions.For numerical reasons, we define r 0 to be at the first outer shell, which has a concentration ∑ T c T (r 0 ) > 0.1, and set the oxygen pressure at ρ(r 0 ) = ρ 0 = 100 mmHg.This assumption produces relatively similar results for the necrotic radius r an when compared to the equation of Grimes et al. [24] if the tumor volume is estimated by r ∑ T c T r i V i (Figure A9).For simplicity, we evaluate the step function at the right edge of each shell i: r * i = (i + 1) • dr.This enables us to estimate the integrals as sum over each shell at the discrete points r * i .We used the following discretization for the integrals of form q 0 x k • c * (x)dx in the above solution of Equation (A18): using the single cell diameter dr * as spatial discretization, which is smaller than the shell width dr and as such increases the spatial resolution of the oxygen profile.While it is in principle possible to set the resolution even smaller than a single cell diameter dr * , our fit results (Table 1) suggest that a spatial resolution smaller than a single cell diameter is not relevant for experimental setups.We translate the cell concentration c T (r i ) to this finer discretization c * (r) via linear interpolation.Due to our missing third boundary condition, ∂ρ(r an ) ∂r = 0, we do not yet know the exact value of the anoxic radius r an ; therefore, we implement a hybrid-analytical approach by using a bisection algorithm to iterate r an over the range [0, r 0 ] until the last boundary condition is met.In detail, we found that if r an is set too high then the minima for ρ * (r) is below r an while if r an is set too low then the minima of ρ * (r) is higher than r an (Figure A10), which we used for the bisection algorithm. .The oxygen profile estimated with the hybrid analytical approach (solid blue) is similar to the previously published method of Grimes et al. [24] (solid black), which assumes hard spheres for the spheroid and its necrotic core.Shown here is a cell concentration for proliferation-competent cells (solid green), membrane-defect cells (solid red), and their sum (dotted black), with the parameter fitted for the FaDu cell line.
. Using bisection, the anoxic radius r an for the oxygen profile is determined by requiring the minima of the oxygen profile ρ * (r) to be at r an , with ρ * (r an ) = ρ an .(a) Example cell concentration of oxygen consuming proliferation-competent cells c p (solid dark green) with corresponding oxygen profiles for different anoxic radii r an (dashed lines).We used the parameters calibrated for FaDu.(b) Magnification of the black box in (a) around the minima of the oxygen profiles.The solid vertical lines in corresponding colors highlight the anoxic radius r an for each oxygen profile.The correct radius determined by bisection is shown in red.
Finally, we can assemble the actual oxygen profile ρ for further processing, which reads as follows: ρ an, for 0 ≤ r ≤ r an ρ * i , for r an ≤ r ≤ r 0 ρ 0 , for r 0 ≤ r ≤ r max .

(A20)
For the further simulation of the RS model, we only need the radii r an and r h .The second can be estimated through piecewise linear interpolation of ρ * i .For computation of the oxygen enhancement ratio OER(ρ(r i )) in the LQ model, we need to perform a further extrapolation from ρ * at points r * i to r i .It is known that the oxygen pressure ρ 0 at the surface of the spheroid is not actually constant, and varies with size of the spheroid typically within the range ρ 0 ∈ [80, 120] mmHg [45,46].However, we confirmed that the goodness of fit of our model prediction is insensitive to variations of ρ 0 within this range (Tables A1 and A2).
Table A1.The RS model is sufficiently insensitive to the oxygen level ρ 0 assumed at the surface of the spheroid (radius r 0 ).Here, we show that for the range ρ 0 ∈ [80, 120] mmHg, which can be expected in experiments [45,46], the goodness of fit is only slightly affected.Examples show the cell lines HCT-116 and MDA-MB-468; the first row shows the fit results of Table 1.Below are two rows where we respectively adjusted ρ 0 and kept the previous fitted values, while separated by the dashed line are two more rows in which we respectively adjusted ρ 0 and refitted the parameters.

Cell Line
HCT-116 Figure 2 100 (main text) Table A2.The RS model is sufficiently insensitive to the oxygen level ρ 0 assumed at the surface of the spheroid (radius r 0 ).Here, we show two that for the range ρ 0 ∈ [80, 120] mmHg, which can be expected in experiments [45,46], the goodness of fit is only slightly affected.Examples show the fit results for radiation doses of 5 Gy and 10 Gy for the cell line HCT-116, where the first row shows the fit results of Table 2. Below are two rows where we respectively adjusted ρ 0 and kept the previous fitted values while separated by the dashed line two more rows where we respectively adjusted ρ 0 and refitted P 1 mc and P 2 mc .All other parameters are taken from calibration of untreated growth (Table 1).When all cells of a spheroid are damaged after treatment, the spheroid dynamics are expressed by The model always exhibits exponential behavior for R < dr, as only one shell is occupied.We estimate an upper bound beyond which the model is fully in the linear regime, assuming that the transition from exponential to linear is smooth without any local maxima: By substituting this time point into Equation (A29), the radius at which the exponential regime has the same first derivative as the linear regime is obtained as R spheroid = 3 • dr.Thus, the transition between both regimes is expected to be in the range dr ≤ R spheroid ≤ 3 • dr.
We can make a further prediction for the necrotic radius if no more proliferationcompetent cells become membrane-defect cells and if the volume occupied by membranedefect cells is reduced with rate δ.In this case, the same Equations (A24)-(A26) can be applied, with the rate γ replaced by δ, to describe the dynamics of membrane-defect cells instead of proliferation-competent cells A comparison of these analytical predictions and our simulation can be found in Figure A11, where we used the same parameters fitted for the cell line FaDu and a radiation dose of 20 Gy while slightly modifying the experimental setup to emphasize the characteristics of the spheroid dynamics.The outer radius R spheroid of the RS model agrees with the analytical predictions.We used the parameter fitted for the FaDu cell line and a radiation dose of 20 Gy.Certain parameters were altered to make the characteristics of the spheroid dynamics more pronounced: the initial volume is lower and the probability for a mitotic catastrophe P mc changes after t = 25 d from P 1 mc = 0 to P 2 mc = 0.87.From the analytical derivation, a linear scaling (black dashed lines) for larger and an exponential scaling (black dotted lines) for smaller radii is expected; (a,b) display the same results in the linear and semi-log-y scales, respectively.The gray area marks the predicted transition region between the linear and exponential regimes dr ≤ R spheroid ≤ 3 • dr.
Note that our analysis shows only an upper bound for the actual growth behavior of the model, which can be seen in Figure A11.There is a smooth transition of the expo-nential growth to the linear depending on how much larger the spheroid radius R spheroid is compared to the shell width dr.This implies that the radial growth is only fully linear for R spheroid → ∞ and fully exponential for R spheroid → 0, resulting in transient behavior between these two extremes.Our findings are in line with previous reported results [65].In conclusion, it is possible to fit various pairs of proliferation rates γ and shell widths dr to reproduce the spheroid growth in the linear regime, which is the one usually explored in experiments.However, if spheroids are observed at a smaller radius R spheroid < 3dr, the transition to the exponential regime should be visible, allowing γ and dr to be uniquely identified.We demonstrate the sensitivity to these two parameters through an example for the FaDu cell line with 20 Gy irradiation in Figure A12, where the dynamics resulting from three different choices for γ and dr can only be discriminated below R spheroid < 150 µm.The dynamics resulting from different choices of γ and dr can only be discriminated at small R spheroid < 150 µm, where the transition between the linear and exponential regimes is observable.As an example, three choices of parameter sets with the same product γdr (same slope in the linear regime) are displayed for the FaDu cell line with 20 Gy irradiation (the solid line corresponds to parameters in Figure A6).The time ln(2)/γ is varied in the fitting range used for the calibration of the RS model (20 h dashed line and 40 h dotted line).

Appendix G. Computational Speed Comparison
The newly formulated RS model can simulate the growth of a tumor spheroid much faster than cell-based models.This fact is highlighted in Figure A13, where we compare the computational time needed for an equal set of parameters to reach a certain tumor spheroid volume for a cell-based model, the new RS model, and a nonspatial model, here, the model of Grimes et al. [24].This leads to the consequence that most of the time cells are not able to actually move inwards due to missing free space when the inward transport rate is higher than the speed at which empty space arises.This results in computational effort without making any simulation progress.To avoid this in the cellular automaton.we implement the concept of a rate-independent inwards shuffle; after every MCS, each cell is allowed to move inwards if free space exist within their interaction neighborhood.Cells located further inwards are allowed to migrate primarily to cells that are further away from the center.When following this algorithm, the tumor spheroid in the cellular automaton stays relatively dense and the event of inwards transportation can be neglected.This results in the following revised equation for the real time step length dt: Note that due to the implementation of the inwards shuffle the parameter mapping from the RS model to this cellular automaton is only successful if the tumor spheroid in both models is relatively dense.In the case where the spheroid in the cellular automaton is not dense, it is necessary to decrease dt in Equation (A37) e.g., by increasing the factor 10.However, this was not necessary in the example shown here.
Because there is no inwards transport in the cellular automaton anymore, the remaining events for each cell type are never simultaneously possible due to the necessary prerequisites, i.e., a proliferation-competent cell T = p can either proliferate if the oxygen level ρ(x, t) is sufficiently large, or die and become anoxic if the oxygen profile is to low; it can never be in a position where it is able to both proliferate and die.Therefore, we have never to choose between two events when we know the cell type and the oxygen level at its position.Thus, each drawn cell within an MCS undergoes the following algorithm: 1.
For each possible event for the selected cell type T, find the one where the prerequisites are fulfilled given the oxygen level ρ(x, t) at the position of the cell and the free space in its neighborhood 2.
If a suitable event is found, When the rate r ate of that specific event follows: r ate • dt > r andom , perform the event.

3.
Otherwise, the drawn cell does not perform any event.
Additionally, we implement an algorithm for calculating the 3D oxygen profile ρ(x) corresponding to the new cell configuration η t+1 within the tumor spheroid after each MSC.Two important things need to be updated.(1) As the tumor grows, we need to re-evaluate the outer surface of the spheroid where the radius r 0 is located and at which the oxygen level ρ 0 is defined, meaning that we keep the oxygen level at the surface of the spheroid constant during the simulation, as in the RS model.(2) We need to update the oxygen profile within in tumor spheroid based on the new radius r 0 and the new specific cell configuration η t+1 .
In step (1), we calculate the size of the spheroid and set the oxygen pressure ρ 0 at the surface of the spheroid.For comparability reasons, we use the same algorithm as in the RS model.Therefore, we divide the 3D tumor into spherical shells of width dr, the same value as in the RS model, and for each shell we calculate the cell concentration per type T ∈ {p, n}.This results in a 1D cell concentration profile similar to the one in the RS model.Next, we can apply the same algorithm used in the RS model to estimate r 0 ; for more details, see Appendix D.
In step (2), we calculate the 3D oxygen profile within the spheroid using a numerical approach with discrete time steps that stop if a steady state is reached.Because the time scale for oxygen diffusion (∼ seconds) is orders of magnitude smaller than the time scale of cell cycle progression (∼ hours), it is reasonable that the steady state is reached on a time-scale shorter than an MCS.The numerical approach is based on Equations (A11) and (A12), and Table A3.The fitted shell size κ of the RS model can be mapped to different neighborhood environments of a cell-based model when transferring the parameter to achieve a similar simulation.If two neighborhood condition are shown, e.g., 1-von Neumann/1-Moore, one is chosen randomly for each successive proliferation event.Note that the 3D cellular automaton itself can have variances for the radius R spheroid that influence the κ fitting.Therefore, we used the mean of five 3D cellular automaton simulations to reduce this error.The correlation is independent of the chosen proliferation rate γ for the range of experimentally expected values.Note that the variance for each fit of κ was smaller than the shown symbols.(b) It is possible to find a linear correlation between the shell size κ and the mean distance at which a cell can be placed from its mother cell in a certain neighborhood condition in the 3D cellular automaton.Note that neum1moor1 refers to a random choice between these two neighborhoods for each proliferation event.The mean radius of R spheroid of 5 simulations of the 3D cellular automaton where used to fit the κ of the RS model.The fitted parameters of the FaDu cell line were used.Only and δ were set to 0 in order to single out the effect of proliferation.Additionally, the starting volume was multiplied by 2.5 to increase visibility.
To compare simulations with experimental data, similar starting conditions are important.However, for the experiments presented here, only the starting volume and respective radii Rspheroid are known, not the radial cell distribution of proliferation-competent or membrane-defect cells.Therefore, we initialized our model with a sigmoid function accounting for a share θ of the starting volume.Then, we simulated until reaching 100% of the starting volume and let the dynamics of the model determine the cell distributions, and consequently the cell concentrations which best agree with the starting radius.
For small θ, e.g., 20%, this is quite simple, as in the realm of the experiments used here only proliferation-competent cells are present.However, to improve the performance of the parameter fitting, we set θ to a higher value, e.g., 90%.This results in much shorter simulation times being required to reach the desired 100%, as the starting condition needs to be calculated for each parameter set independently.However, when starting with a volume of θ = 90%, in certain cases there were membrane-defect cells present, which increases complexity, as the exact structures of the proliferation-competent and membrane-defect cell concentrations are unknown.
Therefore, we extended the above algorithm slightly.In detail, we pre-described a steep sigmoid function and shifted the function along the radius with step size dr.For each position, we then calculated the resulting volume, finally choosing the position which yielded the closest volume to θ = 90% of V 0 .This function returns a cell concentration c that does not differentiate between cell types T. Next" following Grimes et al. [24], we estimated the volume V n of membrane-defect cells for a spheroid with θ = 90% of the initial volume V 0 using Equation (A42) below: and applied the same algorithm as above.This yields the radial profile of membranedefect cells c n and the concentration of proliferation-competent cells equaling c p = c − c n .Subsequently, we simulated the tumor growth until it reached 100% of the starting volume V 0 .During this transient, the shape of the logistic function can relax to its appropriate form, then the given cell concentration at this point can be set as the initial condition for the membrane-defect and proliferation-competent cells.Figure A15 shows an example verifying that our results are independent of the fraction θ ∈ {20%, 90%} of the initial volume used to initialize the simulation. .The simulation result is independent of the starting condition.We compared two different starting cell concentrations for the radial-shell model.The brighter colors (solid) are from the simulation where the starting condition was set as described in Appendix J and used in all other simulations of the manuscript.The darker colored lines (dashed) are from a simulation in which the initial cell concentration was set with a relatively steep sigmoid function for 20% of the starting volume.The same parameters were used for both simulations.

Figure 1 .
Figure 1.Illustration of the radial-shell model (RS-model).(a)Sketch of a typical stained histological section of a spheroid with a viable rim (green) surrounding a secondary necrotic core, with cell debris membrane-defect cells (red) suggesting an approximate rotational symmetry of the spheroid.(b) By exploiting the symmetry, the model is able to consider the cell dynamics on radial shells i ∈ [0, 1, 2, ...) and along the distance r from the center of the spheroid at i = 0.The shells have equal radial width dr, and their maximal volume V i is consequently increasing outwards.Each shell i holds a concentration of cells c T (r i ) of type T cells, where T = p denotes proliferation-competent cells (green, proliferate with maximal rate γ, consume oxygen with rate a) and T = n membranedefect cells (red, irreversibly lost competence to proliferate, do not consume oxygen).Note that the concentrations c T ∈ [0, 1] are normalized according to the available space in each shell.(c) Oxygen consumption of proliferation-competent cells leads to a radial decrease of oxygen pressure from the value ρ 0 at the most outer shell at r 0 down to an anoxic threshold ρ an at r an towards the center of the spheroid.Beyond the threshold ρ an , cells are anoxic and become membrane-defect with the anoxic death rate .The volume occupied by membrane-defect cells is reduced with rate δ and all cells are transported inward with rate λ such that the spheroid remains compact.Then, the outer radius R spheroid and necrotic radius R necrotic are defined as the radius of a sphere with a volume equal to the total cell volume and membrane-defect cell volume, respectively.For reference, the shell width dr = κdr * is expressed as multiple κ of a single cell diameter dr * .

Figure 2 .
Figure 2. (a)The RS model reproduces the outer radius (dark green) for untreated growth of the cell line HCT-116, with experimental data from Brüningk et al.[1] (black crosses) even better than the model proposed along with the data (light green dashed).Note that the necrotic radius during the experiment can be estimated (gray crosses; see Section 3 for details) and that the RS model additionally reproduces this necrotic radius (red).(c) Volumes corresponding to the radii displayed in (a) in a semi-log plot.Multiples of 2 and 5 of the starting volume V t=0 are highlighted with horizontal gray dashed and gray dotted lines, respectively.The RS model predicts the cell concentration over the radius at any time point.Examples show the cell distributions for (b) the first time point t = 0 d and (d) the last day t = 21 d of the experiment.Displayed are the cell concentrations of proliferationcompetent cells c p (green solid), membrane-defect cells c n (red solid), total concentration of cells (black dotted), and oxygen profile (blue solid, right y axis).The marker at the bottom of the panels (b,c) indicates the necrotic and outer radius estimated from the model (red and green crosses) and from the actual experimental data (gray and black arrows), respectively.The used parameters are reported in Table1.
exemplifies the influence of different pairs of the proliferation rate γ and parameter κ for the FaDu cell line with a 20 Gy irradiation dose.

Figure 3 .
Figure 3. (a)The RS model reproduces the outer radius (dark green) for the cell line HCT-116 with 10 Gy radiation experimental data from Brüningk et al.[1] (black crosses) as effectively as the model proposed along with the data (light green dashed).We used the parameters calibrated for untreated growth, only fitted the probabilities of mitotic catastrophe P 1 mc , P 2 mc , and adopted the radiosensitivities α RT and β RT from Brüningk et al.[1].Note that for this regime no secondary necrotic core is expected (experimental estimate shown as gray crosses), which is reproduced by the RS model (red).The parameters we used can be found in Tables1 and 2. (c) Volumes corresponding to the radii displayed in (a) in a semi-log plot.Multiples of 2 and 5 of the starting volume V t=0 are highlighted with horizontal gray dashed and gray dotted lines, respectively.The RS model predicts the cell concentration over the radius at any time point.Examples show the cell distribution for (b) the first time point t = 0 d and (d) the last day t = 21 d of the experiment.Displayed are the cell concentrations of proliferation-competent cells c p (green solid) and membrane-defect cells c n (red solid), the total concentration of cells (black dotted), and the oxygen profile (blue solid, right y axis).The marker at the bottom of the panels (b,c) indicates the necrotic and outer radius estimated from the model (red and green crosses) and from the actual experimental data (gray and black arrows), respectively.

Figure 4 .
Figure 3. (a)The RS model reproduces the outer radius (dark green) for the cell line HCT-116 with 10 Gy radiation experimental data from Brüningk et al.[1] (black crosses) as effectively as the model proposed along with the data (light green dashed).We used the parameters calibrated for untreated growth, only fitted the probabilities of mitotic catastrophe P 1 mc , P 2 mc , and adopted the radiosensitivities α RT and β RT from Brüningk et al.[1].Note that for this regime no secondary necrotic core is expected (experimental estimate shown as gray crosses), which is reproduced by the RS model (red).The parameters we used can be found in Tables1 and 2. (c) Volumes corresponding to the radii displayed in (a) in a semi-log plot.Multiples of 2 and 5 of the starting volume V t=0 are highlighted with horizontal gray dashed and gray dotted lines, respectively.The RS model predicts the cell concentration over the radius at any time point.Examples show the cell distribution for (b) the first time point t = 0 d and (d) the last day t = 21 d of the experiment.Displayed are the cell concentrations of proliferation-competent cells c p (green solid) and membrane-defect cells c n (red solid), the total concentration of cells (black dotted), and the oxygen profile (blue solid, right y axis).The marker at the bottom of the panels (b,c) indicates the necrotic and outer radius estimated from the model (red and green crosses) and from the actual experimental data (gray and black arrows), respectively.

Figure 5 .
Figure 4. (a)The RS model reproduces the outer radius (dark green) for the cell line HCT-116, with 5 Gy radiation experimental data from Brüningk et al.[1] (black crosses) even better than the model proposed along with the data (light green dashed).We used the parameters calibrated for untreated growth, the values for P 1 mc and P 2 mc fitted at 10 Gy, and the values α RT and β RT from Brüningk et al.[1].Note that the necrotic radius during the experiment can be estimated (gray crosses) (see Section 3 for details) and that the RS model reproduces this necrotic radius (red).Additionally, note that the necrotic radius of the RS model appears jagged due to the spatial discretization of r 0 .The used parameters are reported in Tables1 and 2. (c) Volumes corresponding to the radii displayed in (a) in a semi-log plot.Multiples of 2 and 5 of the starting volume V t=0 are highlighted with horizontal gray dashed and gray dotted lines, respectively.The RS model predicts the cell concentration over the radius at any time point.Examples show the cell distribution for (b) the first time point t = 0 d and (d) the last day t = 21 d of the experiment.Displayed are the cell concentrations of proliferationcompetent cells c p (green solid) and membrane-defect cells c n (red solid), the total concentration of cells (black dotted), and the oxygen profile (blue solid, right y axis).The marker at the bottom of the panels (b,c) indicates the necrotic and outer radius estimated from the model (red and green crosses) and from the actual experimental data (gray and black arrows), respectively.

Figure 6 .
Figure6.The cellular automaton reproduces the experimental growth curve for the FaDu cell line after parameter transfer from the previously fitted RS model.We directly transferred the parameters of the RS model to the cellular automaton using a 50:50 alteration between the 3-Moore and 4-Moore neighborhoods to achieve a value of κ = 4.26.(a) Cell concentration over the radius at t = 14 d for the 3D cellular automaton; the mean cell concentrations (solid) from ten independent simulation runs are compared with the cell concentrations of the corresponding RS model (dashed).The cell concentrations of proliferation-competent cells c p (green solid) and membrane-defect cells c n (red solid) are displayed.(b) Growth dynamics over time for the outer radius R spheroid (green) and necrotic radius R necrotic (red) along with the experimental outer radius (black crosses) and corresponding estimated necrotic radius (gray crosses).For the cellular automaton, the maximal and minimal radii from ten independent simulations are shown.
Figure A1.(a)The RS model reproduces the outer radius (dark green) for untreated growth of the cell line HCT-116, with experimental data from Grimes et al.[24] (black crosses) even better than the model proposed along with the data (light green dashed).Note that the necrotic radius during the experiment can be estimated (gray crosses; see main text Section 3 for details) and that the RS model additionally reproduces this necrotic radius (red).(c) Volumes corresponding to the radii displayed in (a) in a semi-log plot.Multiples of 2 and 5 of the starting volume V t=0 are highlighted with horizontal gray dashed and gray dotted lines, respectively.Additionally, the RS model predicts the cell concentration over the radius at any time point.Examples show the cell distribution for (b) the first time point t = 0 d and (d) the last day t = 17 d of the experiment.Displayed are the cell concentrations of proliferation-competent cells c p (green solid) and membrane-defect cells c n (red solid), the total concentration of cells (black dotted), and the oxygen profile (blue solid, right y axis).The marker at the bottom of the panels (b,c) indicates the necrotic and outer radius estimated from the model (red and green crosses) and from the actual experimental data (gray and black arrows), respectively.The used parameters are reported in Table1.

Figure A2 .
Figure A2.(a)The RS model reproduces the outer radius (dark green) for untreated growth of the cell line MDA-MB-468, with experimental data from Grimes et al.[24] (black crosses) even better than the model proposed along with the data (light green dashed).Note that the necrotic radius during the experiment can be estimated (gray crosses; see main text Section 3 for details) and that the RS model additionally reproduces this necrotic radius (red).(c) Volumes corresponding to the radii displayed in (a) in a semi-log plot.Multiples of 2 and 5 of the starting volume V t=0 are highlighted with horizontal gray dashed and gray dotted lines, respectively.Additionally, the RS model predicts the cell concentration over the radius at any time point.Examples show the cell distribution for (b) the first time point t = 0 d and (d) the last day t = 12 d of the experiment.Displayed are the cell concentrations of proliferation-competent cells c p (green solid) and membrane-defect cells c n (red solid), the total concentration of cells (black dotted), and the oxygen profile (blue solid, right y axis).The marker at the bottom of the panels (b,c) indicates the necrotic and outer radius estimated from the model (red and green crosses) and from the actual experimental data (gray and black arrows), respectively.The used parameters are reported in Table1.

Figure A3 .
Figure A3.The RS model reproduces the outer radius (dark green) for untreated growth of the cell line LS-174T, with experimental data from Grimes et al.[24] (black crosses) even better than the model proposed along with the data (light green dashed).Note that the necrotic radius during the experiment can be estimated (gray crosses; see main text Section 3 for details) and that the RS model additionally reproduces this necrotic radius (red).(c) Volumes corresponding to the radii displayed in (a) in a semi-log plot.Multiples of 2 and 5 of the starting volume V t=0 are highlighted with horizontal gray dashed and gray dotted lines, respectively.Additionally, the RS model predicts the cell concentration over the radius at any time point.Examples show the cell distribution for (b) the first time point t = 0 d and (d) the last day t = 7 d of the experiment.Displayed are the cell concentrations of proliferation-competent cells c p (green solid) ans membrane-defect cells c n (red solid), the total concentration of cells (black dotted), and the oxygen profile (blue solid, right y axis).The marker at the bottom of the panels (b,c) indicates the necrotic and outer radius estimated from the model (red and green crosses) and from the actual experimental data (gray and black arrows), respectively.The used parameters are reported in Table1.

25 Figure A4 .
Figure A4.The RS model reproduces the outer radius (dark green) for untreated growth of the cell line SCC-25, with experimental data from Grimes et al.[24] (black crosses) even better than the model proposed along with the data (light green dashed).Note that the necrotic radius during the experiment can be estimated (gray crosses; see main text Section 3 for details) and that the RS model additionally reproduces this necrotic radius (red).(c) Volumes corresponding to the radii displayed in (a) in a semi-log plot.Multiples of 2 and 5 of the starting volume V t=0 are highlighted with horizontal gray dashed and gray dotted lines, respectively.Additionally, the RS model predicts the cell concentration over the radius at any time point.Examples show the cell distribution for (b) the first time point t = 0 d and (d) the last day t = 17 d of the experiment.Displayed are the cell concentrations of proliferation-competent cells c p (green solid) and membrane-defect cells c n (red solid), the total concentration of cells (black dotted), and the oxygen profile (blue solid, right y axis).The marker at the bottom of the panels (b,c) indicates the necrotic and outer radius estimated from the model (red and green crosses) and from the actual experimental data (gray and black arrows), respectively.The used parameters are reported in Table1.

FaDuFigure A5 .
Figure A5.The RS model reproduces the outer radius (dark green) for untreated growth of the FaDu cell line, with experimental data from Chen et al.[18] (black crosses).Note that the necrotic radius during the experiment can be estimated (gray crosses; see main text Section 3 for details) and that the RS model additionally reproduces this necrotic radius (red).(c) Volumes corresponding to the radii displayed in (a) in a semi-log plot.Multiples of 2 and 5 of the starting volume V t=0 are highlighted with horizontal gray dashed and gray dotted lines, respectively.Additionally, the RS model predicts the cell concentration over the radius at any time point.Examples show the cell distribution for (b) the first time point t = 0 d and (d) the last day t = 14 d of the experiment.Displayed are the cell concentrations of proliferation-competent cells c p (green solid) and membrane-defect cells c n (red solid), the total concentration of cells (black dotted), and the oxygen profile (blue solid, right y axis).The marker at the bottom of the panels (b,c) indicates the necrotic and outer radius estimated from the model (red and green crosses) and from the actual experimental data (gray and black arrows), respectively.The used parameters are reported in Table1.

FaDuFigure A6 .
Figure A6.(a)The RS model reproduces the outer radius (dark green) for the FaDu cell line with 20 Gy radiation, experimental data from Chen et al.[18] (black crosses) as effectively as the model proposed along with the data (light green dashed).We used the parameters calibrated for untreated growth, only fitted the probabilities of mitotic catastrophe P 1 mc , P 2 mc , and adopted the radiosensitivities α RT and β RT from the ranges observed in Chen et al.[17] and Xu et al.[59].Note that no necrotic core is expected for this regime (experimental estimate shown as gray crosses), which is reproduced by the RS model (red).The used parameters can be found in Tables1 and 2. (c) Volumes corresponding to the radii are displayed in (a) in a semi-log plot.Multiples of 2 and 5 of the starting volume V t=0 are highlighted with horizontal gray dashed and gray dotted lines, respectively.Additionally, the RS model predicts the cell concentration over the radius at any time point.Examples show the cell distribution for (b) the first time point t = 0 d and (d) the last day t = 33 d of the experiment.Displayed are the cell concentrations of proliferation competent cells c p (green solid) and membranedefect cells c n (red solid), the total concentration of cells (black dotted), and the oxygen profile (blue solid, right y axis).The marker at the bottom of the panels (b,c) indicates the necrotic and outer radius estimated from the model (red and green crosses) and from the actual experimental data (gray and black arrows), respectively.
Figure A7.(a)The RS model reproduces the outer radius (dark green) for the FaDu cell line with 5 Gy radiation, experimental data from Chen et al.[18] (black crosses).We used the parameters calibrated for untreated growth, the values for P 1 mc and P 2 mc fitted at 20 Gy, and adopted the radiosensitivities α RT and β RT from the ranges observed in Chen et al.[17] and Xu et al.[59].Note that no necrotic core is expected for this regime (experimental estimate shown as gray crosses), which is reproduced by the RS model (red).The used parameters can be found in Tables1 and 2. (c) Volumes corresponding to the radii are displayed in (a) in a semi-log plot.Multiples of 2 and 5 of the starting volume V t=0 are highlighted with horizontal gray dashed and gray dotted lines, respectively.Additionally, the RS model predicts the cell concentration over the radius at any time point.Examples show the cell distribution for (b) the first time point t = 0 d and (d) the last day t = 21 d of the experiment.Displayed are the cell concentrations of proliferation competent cells c p (green solid) and membranedefect cells c n (red solid), the total concentration of cells (black dotted), and the oxygen profile (blue solid, right y axis).The marker at the bottom of the panels (b,c) indicates the necrotic and outer radius estimated from the model (red and green crosses) and from the actual experimental data (gray and black arrows), respectively.
Figure A8.(a)The RS model reproduces the outer radius (dark green) for the FaDu cell line, with 2.5 Gy radiation experimental data from Chen et al.[18] (black crosses).We used the parameters calibrated for untreated growth, the values P 1 mc and P 2 mc fitted at 20 Gy, and adopted the radiosensitivities α RT and β RT from the ranges observed in Chen et al.[17] and Xu et al.[59].The necrotic radius of the experiments can be estimated (gray crosses; see Section 3 for details).The RS model can predict the necrotic radius (red) in a good manner.Note that the necrotic radius of the RS model appears jagged due to the spatial discretization of r 0 .The used parameters can be found in Tables1 and 2. (c) Shows the same data plotted for the volume in a semi-log plot.Multiples of 2 and 5 of the starting volume V t=0 are highlighted with horizontal gray dashed and gray dotted lines, respectively.Additionally, the RS model predicts the cell concentration over the radius at any time point.Examples show the cell distribution for (b) the first time point t = 0 d and (d) the last day t = 14 d of the experiment.Displayed are the cell concentrations of proliferation competent cells c p (green solid) and membrane-defect cells c n (red solid), the total concentration of cells (black dotted), and the oxygen profile (blue solid, right y axis).The marker at the bottom of the panels (b,c) indicates the necrotic and outer radius estimated from the model (red and green crosses) and from the actual experimental data (gray and black arrows), respectively.
Figure A9.The oxygen profile estimated with the hybrid analytical approach (solid blue) is similar to the previously published method of Grimes et al.[24] (solid black), which assumes hard spheres for the spheroid and its necrotic core.Shown here is a cell concentration for proliferation-competent cells (solid green), membrane-defect cells (solid red), and their sum (dotted black), with the parameter fitted for the FaDu cell line.
Figure A11.The outer radius R spheroid of the RS model agrees with the analytical predictions.We used the parameter fitted for the FaDu cell line and a radiation dose of 20 Gy.Certain parameters were altered to make the characteristics of the spheroid dynamics more pronounced: the initial volume is lower and the probability for a mitotic catastrophe P mc changes after t = 25 d from P 1 mc = 0 to P 2 mc = 0.87.From the analytical derivation, a linear scaling (black dashed lines) for larger and an exponential scaling (black dotted lines) for smaller radii is expected; (a,b) display the same results in the linear and semi-log-y scales, respectively.The gray area marks the predicted transition region between the linear and exponential regimes dr ≤ R spheroid ≤ 3 • dr.

Figure A12.
Figure A12.The dynamics resulting from different choices of γ and dr can only be discriminated at small R spheroid < 150 µm, where the transition between the linear and exponential regimes is observable.As an example, three choices of parameter sets with the same product γdr (same slope in the linear regime) are displayed for the FaDu cell line with 20 Gy irradiation (the solid line corresponds to parameters in FigureA6).The time ln(2)/γ is varied in the fitting range used for the calibration of the RS model (20 h dashed line and 40 h dotted line).

Figure A14 .
Figure A14.The shell size κ of the RS model can be related to the individual neighborhood condition of a 3D cellular automaton.(a) The correlation is independent of the chosen proliferation rate γ for the range of experimentally expected values.Note that the variance for each fit of κ was smaller than the shown symbols.(b) It is possible to find a linear correlation between the shell size κ and the mean distance at which a cell can be placed from its mother cell in a certain neighborhood condition in the 3D cellular automaton.Note that neum1moor1 refers to a random choice between these two neighborhoods for each proliferation event.The mean radius of R spheroid of 5 simulations of the 3D cellular automaton where used to fit the κ of the RS model.The fitted parameters of the FaDu cell line were used.Only and δ were set to 0 in order to single out the effect of proliferation.Additionally, the starting volume was multiplied by 2.5 to increase visibility.
Figure A15.The simulation result is independent of the starting condition.We compared two different starting cell concentrations for the radial-shell model.The brighter colors (solid) are from the simulation where the starting condition was set as described in Appendix J and used in all other simulations of the manuscript.The darker colored lines (dashed) are from a simulation in which the initial cell concentration was set with a relatively steep sigmoid function for 20% of the starting volume.The same parameters were used for both simulations.