The evolution of tumour composition during fractionated radiotherapy: implications for outcome

Current protocols for delivering radiotherapy are based primarily on tumour stage and nodal and metastases status, even though it is well known that tumours and their microenvironments are highly heterogeneous. It is well established that the local oxygen tension plays an important role in radiation-induced cell death, with hypoxic tumour regions responding poorly to irradiation. Therefore, to improve radiation response, it is important to understand more fully the spatiotemporal distribution of oxygen within a growing tumour before and during fractionated radiation. To this end, we have extended a spatially-resolved mathematical model of tumour growth first proposed by Greenspan (Stud. Appl. Math., 1972) to investigate the effects of oxygen heterogeneity on radiation-induced cell death. In more detail, cell death due to radiation at each location in the tumour, as determined by the well-known linear-quadratic model, is assumed also to depend on the local oxygen concentration. The oxygen concentration is governed by a reaction-diffusion equation that is coupled to an integro-differential equation that determines the size of the assumed spherically-symmetric tumour. We combine numerical and analytical techniques to investigate radiation response of tumours with different intra-tumoral oxygen distribution profiles. Model simulations reveal a rapid transient increase in hypoxia upon re-growth of the tumour spheroid post-irradiation. We investigate the response to different radiation fractionation schedules and identify a tumour-specific relationship between inter-fraction time and dose per fraction to achieve cure. The rich dynamics exhibited by the model suggest that spatial heterogeneity may be important for predicting tumour response to radiotherapy for clinical applications.


Introduction
Cancer continues to be one of the main causes of mortality worldwide, with 8.2 million cancerrelated deaths estimated to have occurred globally in 2012 (Torre et al, 2015). Of the millions of people diagnosed with some form of cancer each year, about half will receive radiotherapy as part of their treatment (Fowler, 2006).
Typically, total radiation dose is delivered to the tumour as a series of small doses, or fractions, administered over a period of several days or weeks in order to limit the toxic side effects to healthy cell populations. The conventional fractionation schedule for most tumors, where a dose of approximately 2 Gy (Gray) is delivered once a day Monday to Friday up to a total of 50-70 Gy, has remained the standard of care for many years (Ahmed et al, 2014;Marcu, 2010). Whilst altered schemes, such as hyperfractionation, accelerated fractionation and hypofractionation, have been suggested as alternatives for certain indications, the selection of an optimal fractionation protocol for a particular tumour would clearly benefit from a more individualised approach. In particular, beyond tumour location and stage, patient-specific factors that may be important in determining response to a particular treatment are currently not considered. It is in this regard that mathematical modelling has the potential to play an important role; identifying key factors that determine treatment outcomes and ultimately identifying patient-specific treatment protocols.
Radiation damages to a DNA fragment can be either a single-or double-strand break (Hendry, 1979;Thrall, 1997;Muriel, 2006). A normal cell has repair mechanisms designed to correct for single-strand breaks. However these can fail, with misrepair potentially leading to cell death. Whilst double-strand breaks are rarer than single-strand breaks, they are more likely to result in cell death due to increased difficulty in repairing the damaged DNA. Tumours differ in their responses to irradiation, with radiosensitivity being understood as an intrinsic property of the cell population that could be estimated from molecular analysis of biopsy samples (Eschrich et al, 2009). Other factors that affect a tissue's response to radiotherapy include repopulation between radiation doses, tissue re-oxygenation and redistribution of cells within the cell cycle (Thrall, 1997). For rapidly proliferating tumours in particular, repopulation between fractions of radiation is important. For tumours containing larger regions of hypoxia, re-oxygenation of the tumour may play a pivotal role. Ionising radiation interacts with oxygen and water molecules, creating oxygen free radicals that are then able to cause DNA damage (indirect action). Thus, tumour cells in hypoxic regions exhibit decreased response to radiotherapy.
Due to the difficulties of studying tumours in vivo, such effects are often investigated experimentally using in vitro multicellular tumour spheroids. Tumour spheroids provide a controlled study environment of intermediate complexity between 2D culture media and in vivo models in which important insight can be gained into the biology of the tumour micro-environment and the effect of therapies on a growing tumour (Carlson et al, 2006;Hirschhaeuser et al, 2010;Mueller-Klieser, 1987;Sutherland et al, 1981). However, tumour spheroid growth is avascular with a typical diameter of less than 5mm, so while useful for calibrating mathematical models, parameters derived from tumour spheroids may need adjustment before they can be used to simulate in vivo tumours.
Clonogenic assays are widely used to determine clonogenic survival for a particular cell type after acute radiation doses (Joiner and van der Kogel (2009), Chapter 4). Measurements of the clonogenic surviving fraction for a single cell type across a range of doses can be used to generate a survival curve. Such curves are of similar shape across most cell types and give rise to the most widely adopted mathematical model of radiotherapy response, the linear-quadratic (LQ) model ( (Joiner and van der Kogel, 2009), Chapter 4). We note that although this model is empirical, mechanistic models have been proposed to explain it ( (Sachs et al, 1997;Nilsson et al, 1990), (Joiner and van der Kogel, 2009) Chapter 4). In the LQ model, the survival fraction, SF d , of tumour cells after a dose, d Gy, of radiation is given by where α (Gy −1 ) and β (Gy −2 ) are intrinsic radiosensitivity parameters (Fowler, 2006;Sachs et al, 1997;Withers, 1999;O'Rourke et al, 2009). Values for α and β may differ significantly between tumours and patients. The α/β (Gy) ratio is often used to characterize the sensitivity of a particular tissue type to fractionation. Values of α/β can fall as low as 1 Gy for late-responding tissues such as prostate cancer, and reach as high as 20 Gy for early-responding, rapidly proliferating tissues such as head and neck cancer (Withers, 1999). The LQ model is frequently incorporated into more detailed tumour models as an instantaneous effect. Many different mathematical models and methodologies have been used to study tumour growth. Among the simplest of these are ordinary differential equation (ODE) models that aim to qualitatively capture observed growths. Small, initial tumour volumes often exhibit exponential growth dynamics, followed by a deceleration as tumour growth saturates due to microenvironmental effects, such as limited space and nutrient supply, resulting in sigmoidal growth curves (Sachs et al, 2001;McAneney and O'Rourke, 2007). Popular ODE models include the logistic and Gompertz growth models where the volume-saturation limit is represented by the carrying capacity, K.
Mathematical models are often built on the simplifying assumption that a tumour has a spatially homogeneous composition. Phenomenological models contain few parameters and do not account for complex underlying biological interactions. Whilst such models may provide limited biological insight, the ease with which model parameters may be estimated from limited clinical data makes them attractive for making predictions about radiotherapy response and identifying personalised fractionation protocols in the clinic (Prokopiou et al, 2015).
A variety of more complex, spatially-resolved, models have been proposed in order to provide more mechanistic insight, often reflecting the heterogeneous nature of growing tumours. As a tumour grows, it will typically develop hypoxia and necrosis in regions where oxygen or nutrient supply is inadequate. Multicellular, avascular tumour spheroids are used as in vitro models to study, in a controlled environment, effects observed in in vivo tumours (Carlson et al, 2006;Hirschhaeuser et al, 2010;Mueller-Klieser, 1987;Sutherland et al, 1981). In a similar manner, it is natural to develop mathematical models for avascular tumour spheroids before considering the more complex case of vascular tumour growth.
One of the simplest spatially-resolved models of avascular tumour growth was proposed by Greenspan (Greenspan, 1972). In the Greenspan model, the outer tumour radius evolves in response to a single diffusible nutrient species, commonly taken to be oxygen. Internal free boundaries, decomposing the tumour spheroid into a central necrotic core, an outer proliferating rim and an intermediate hypoxic annulus, are determined by contours of the oxygen profile across the tumour.
Under certain simplifying assumptions, Greenspan's model can be reduced to a coupled system of ODEs and algebraic equations. Most other spatially-resolved continuum models are formulated as systems of partial differential equations (PDEs). Such approaches include frameworks from multiphase modelling and morphoelasticity. Multiphase models consider the tumour environment as a mixture of two or more separate fluid cell populations, or phases. Systems of PDEs are obtained by applying mass and momentum balances to each phase and making suitable constitutive assumptions about their properties and interactions (Breward et al, 2003;Byrne et al, 2003;Preziosi and Tosin, 2009). Morphoelastic models of tumour growth provide a theoretical framework within which to study biological tissues for which growth and elasticity are inter-related (Araujo and McElwain, 2004).
While incorporating more biological detail, these more sophisticated approaches typically involve more parameters that may be difficult to estimate in practice. In this paper we present a simple spatially-resolved model for tumour spheroid growth and response to radiotherapy in order to investigate the effects of spatial heterogeneity. In Section 2 we extend Greenspan's original spatial model for tumour spheroid growth (Greenspan, 1972) to include radiation effects. In Section 3.1 we numerically solve the model equations and discuss key features of the resulting dynamics. We use a combination of further numerical simulations and analysis of the model equations in order to explore the tumour dynamics exhibited by the model.

Model development
In this section, we introduce a spatial model of avascular tumour growth originally developed by Greenspan (Greenspan, 1972). We extend Greenspan's model to account for the effects of radiation and arrive at a new model for tumour response to radiotherapy.

Original growth model by Greenspan
Following Greenspan (Greenspan, 1972), we consider the growth of a radially-symmetric, avascular tumour spheroid in response to a single, growth-rate limiting, diffusible nutrient, here oxygen. For simplicity, we assume that growth inhibition occurs due to nutrient deficiency rather than in response to an inhibitory factor. We denote the outer tumour radius by R(t) and let c(r, t) denote the oxygen concentration at a distance 0 ≤ r ≤ R(t) from the tumour centre. We suppose that the oxygen concentration is maintained at a constant level, c ∞ , on r = R(t). Oxygen diffuses on a much shorter timescale than tumour growth, taking ∼ 10 seconds to diffuse 100 µm by comparison with a tumour growth period measured in days (Greenspan, 1972). As such we assume that the oxygen concentration is in a quasi-steady state. Internal free boundaries at r = R H (t) and r = R N (t) are defined implicitly and denote contours on which the oxygen concentration attains the threshold values c H and c N , respectively. These interfaces decompose a well-developed tumour into a central necrotic core (0 < r < R N ) where c ≤ c N , and an outer, oxygen-rich region (R H < r < R) in which c H < c, these regions being separated by a hypoxic annulus (R N < r < R H ) in which c N < c < c H (see Figure 1 for a schematic).
(a) Image of an in vitro tumour spheroid (adapted from (Folkman and Hochberg, 1973)) highlighting the distinct microenvironments.
(b) Schematic of the corresponding model tumour spheroid geometry. Figure 1: Diagrams highlighting the distinct regions within a tumour spheroid and how they are influenced by the oxygen profile. A central necrotic core (0 < r < R N ) is surrounded by a hypoxic annulus (R N < r < R H ) and an outer proliferating rim (R H < r < R). The moving boundaries r = R N (t), R H (t) and R(t) delineate these regions.
Greenspan's original model (Greenspan, 1972) describes how the dependent variables c(r, t), R(t), R H (t) and R N (t) evolve over time and can be written in dimensionless form as follows: where H(.) is the Heaviside function (H(x) = 1 for x > 0, and H(x) = 0 otherwise), and Γ, λ A , λ N , c ∞ , c H & c N are positive constants. Equations (2)-(8) may be reduced to a system involving a single ODE for R(t) and algebraic equations for R H (t) and R N (t). Analysis of the resulting equations and details of the non-dimensionalisation may be found in (Greenspan, 1972;Byrne, 2012), and for further explanation of the underlying model assumptions see Appendix A.1. Figure  2 shows how the size and composition of the tumour spheroid evolve for a typical simulation of the Greenspan model using the dimensionless parameter values in Table 1    thus, a change in the size and composition of the tumour. The efficacy of the radiation and subsequent cell death depends on the local oxygen concentration. As such, we determine the radiation-induced cell death in each tumour region separately. We denote by R ± , R ± H and R ± N the radii immediately before (-) and after (+) radiotherapy in the proliferating, hypoxic and necrotic tumor compartments.

Radiation-induced cell death and the Linear-Quadratic model
We use the linear-quadratic model (Fowler, 2006;Joiner and van der Kogel, 2009;Enderling et al, 2010) to account for cell kill due to radiotherapy in the well-oxygenated outer rim, R H < r < R, so that the volume survival fraction, SF (d), immediately after a dose d of radiation is given by sf normoxic = [volume of normoxic region after radiotherapy] [volume of normoxic region before radiotherapy] for radiosensitivity parameters α and β. Values of α/β vary markedly, with typical values falling in the range of 3-10 Gy. However extremal values of 1 Gy and 20 Gy have been reported for late-responding tissues such as prostate cancer and early-responding, rapidly proliferating tissues such as head and neck cancer, respectively (Withers, 1999). Typical radiosensitivity parameters for rapidly-proliferating, early-responding tumours are shown in Table 2. When exposed to ionising radiation, the potent oxygen free radicals that form in well-oxygenated regions increase the amount of DNA damage by up to a factor of 3 when compared with hypoxic tumour regions (Alper and Howard-Flanders, 1956). Following Carlson et al. (Carlson et al, 2006), we incorporate the oxygen enhancement ratio, OER, to account for this effect. Radiosensitivity parameters are typically quoted for normoxic conditions. As such we take the common approach and use the OER as a constant factor (taken to be equal to 3 in the presented simulations) that reduces the intrinsic radiosensitivity parameters of tumour cells, α and β, in hypoxic regions. This creates a discontinuity in the response to radiotherapy at r = R H . The volume survival fraction immediately after a dose of radiation is delivered to a population of hypoxic cells is given by sf hypoxic = [volume of hypoxic region after radiotherapy] [volume of hypoxic region before radiotherapy] There is no clear consensus in the literature as to how to modify the 'quadratic' component of cell death, however the form used in Equation (10) affords the interpretation of the OER as the multiplying factor for the dose escalation required under hypoxia in order to achieve the same cell kill as under normoxic conditions. We note that, in practice, this effect is likely to depend continuously on the local oxygen concentration. A continuous functional form for the OER (OER = OER(c)) was proposed by Alper & Howard-Flanders (Alper and Howard-Flanders, 1956). For realistic parameter regimes, simulations using the continuous OER presented in (Alper and Howard-Flanders, 1956) and the discrete OER presented here did not significantly vary (results not shown). For these reasons we restrict attention to the discrete OER stated above.
We assume further that dead material within the necrotic core is unaffected by radiation and so we impose that R + N = R − N . Combining these assumptions we deduce that, for a well-developed tumour spheroid ( Figure  2 loss from hypoxic region which, on rearrangement, yields the following expression for R + in terms of R − , R − H , and R ± N = R N : The oxygen concentration profile associated with the new tumour structure can be determined (Equation (2)) and R H defined implicitly as before (see Equation (4)). We note that while not all hypoxic tumour cells will be killed from a single radiotherapy fraction (sf hypoxic > 0), the instantaneous volume loss due to irradiation and the subsequent reoxygenation of the tumour spheroid may result in a post-radiotherapy tumour composition without a hypoxic region.
Two cases can arise when a well-developed, 3-layer tumour is irradiated (see Figure 3): -irradiated tumour has a necrotic core, but no hypoxic annulus. Since the necrotic core is assumed to be unaffected by radiotherapy, these are the only possible options. We have already described how a tumour spheroid of case (i) responds to radiation (Equation (11)). If a tumour with a pre-radiotherapy composition as in case (ii) is irradiated then the corresponding volume change is given by which, on rearrangement, supplies the following expression for R + in terms of R − and R N

Reconciling pre-and post-radiation tumour growth
In the normal growth regime without radiotherapy, R N is defined implicitly by the oxygen concentration (see Equation (5)). When the tumour volume is reduced due to radiation, the oxygen concentration at the centre of the tumour increases and the location of the interface on which c = c N will shift towards the tumour centre or disappear. As the necrotic core is unaffected by radiation, c(R + N ) > c N and the growth of the tumour spheroid immediately post-radiotherapy does not follow the original, control dynamics.
Between fractions of radiotherapy, repopulation of the tumour occurs. The tumour cells proliferate and die as before, however, whilst c(R N ) > c N no new material is added to the necrotic core, and so the existing necrotic material simply decays at the rate λ A + λ N . In this case R N evolves according to where t f raction is the time of the last fraction delivered for which c(R N ) = c N , and R N (t f raction ) is the radius of the necrotic core upon irradiation. As such, prior to treatment with radiotherapy at t = t f raction the tumour composition is consistent with the control, Greenspan dynamics, while Equation (13) governs the evolution of R N post-irradiation.

Summary: statement of full model (dimensionless)
The tumour growth model combined with fractionated radiotherapy (dose d i Gy delivered at times t = t i , i = 1, 2, ...) can be summarised as follows: 1 4π d dt Continuity conditions for c and ∂c ∂r continuous across r = R H , R N are imposed. Radiotherapy, of dose d i applied at times t = t i , effects an instantaneous volume change. The survival fraction of the normoxic tumour cell population is given by whilst the survival fraction in hypoxic regions is given by The necrotic core of the tumour spheroid is unaffected by radiotherapy. Immediately after a dose of radiation, the outer tumour radius is given by NB: We subtly alter the equation for R H from the original Greenspan model so that R H = R N if the threshold for hypoxia, c H is not reached. This does not change the dynamics in the case of untreated tumour growth, but ensures that the ordering of the tumour boundaries is conserved in the post-radiotherapy regime (0 ≤ R N ≤ R H ≤ R).

Investigation of model behaviour
We now consider the effects of various treatment protocols on a given tumour under this model. In Section 3.1 we solve the model numerically and highlight key features of the simulations for different parameter combinations. We explore these features in more detail in the following sections. We analyse the rapid increase in hypoxia during regrowth of the tumour spheroid post-radiotherapy in Section 3.2. In Section 3.3 we study tumour regrowth following a single dose of radiation, and in Section 3.4 we consider the long-term behaviour of a tumour for different fractionation schedules.

Numerical simulation results
For each tumour growth regime, the oxygen profile c(r, t) may be solved analytically. Equations (14)-(23) may then be reduced to an ODE for R and two algebraic equations for R H and R N , as described in Appendix A.1. We solve the resulting system of equations numerically using a finite difference scheme implemented in MATLAB.
In Figure 4a we present results showing a typical tumour response to a conventional fractionation schedule (2 Gy, 5 days / week) simulated for 6 weeks using the parameter values in Tables 1 & 2. Cell death induced by the first fraction of radiation delivered at t = 1 results in the loss of the hypoxic annulus. While c(R N ) > c N , the necrotic core decays exponentially, as defined by Equation (13). Figure 5a shows how the oxygen concentration in the necrotic core evolves throughout the first week of the radiotherapy protocol. We see the mismatch between c(R N ) and c N following delivery of the first fraction (at t = 1). As the tumour grows between fractions, oxygen concentration in the necrotic core decreases, while each fraction of radiation and the corresponding volume loss results in reoxygenation.
For this parameter combination, radiotherapy results in a tumour volume at the end of treatment that is sufficiently small such that the spheroid is almost entirely composed of proliferative cells. In particular, the necrotic core has decayed so that R N R. Figure 4b shows how the system dynamics change when a tumour with a much lower apoptotic rate (λ A ) is exposed to the standard fractionation protocol. The slower rate of apoptosis results in a larger steady state tumour in the absence of radiotherapy. Similarly, a larger tumour volume supported throughout treatment and a more gradual volume reduction. In this case the hypoxic annulus reappears during treatment (as observed during the weekend breaks in the protocol) since the necrotic core decays such that c(R N ) = c H before the end of the fractionation protocol. A rapid transient increase in R H is observed when the hypoxic annulus reappears. We characterize this behaviour in Section 3.2.
For the simulation in Figure 4b, a tumour spheroid comprising a necrotic core and transient hypoxia persists at the end of treatment. The tumour volume appears to be evolving towards a periodic orbit and, as such, the model predicts that continuing the same fractionation schedule indefinitely will not yield significant further volume reduction.
Since the tumour composition changes dynamically throughout treatment, the efficacy of each radiation fraction also varies. The survival fraction changes by about 10% from the start to the end of treatment due to the shrinkage of the necrotic core (Figure 4b). In this case a characteristic shape within the curve in Figure 5b repeated weekly due to the weekly oscillations in tumour composition, with the spike associated with the rapid reappearance of hypoxia. Such dynamic behaviour clearly depends on tumour-specific parameters and highlights the additional details that can be observed when spatial effects are included in a mathematical model (c.f. difference in tumour composition between Figures 4a & 4b).  Tables 1 & 2 with an initial tumour radius of 3. In (b) we simulate irradiating a tumour spheroid with lower rates of apoptosis (λ A = 0.1213) and oxygen consumption (Γ = 0.3032) starting from an overall radius of 7.5. The evolution of the tumour contours R, R H and R N is shown by the blue, red and yellow lines, respectively.

Asymptotic analysis of model behaviour post-radiotherapy
When a fully-developed tumour spheroid including a hypoxic annulus and necrotic core is exposed to fractionated radiation, the composition of the spheroid immediately after irradiation depends on the dose d, as described in Section 2.2.1. If the dose is sufficiently large, the post-radiotherapy composition comprises a proliferating rim and a necrotic core, but with the absence of a region of hypoxia. In this scenario, the oxygen concentration in the necrotic core is greater than the threshold for hypoxia (i.e. c(R N ) > c H ). Upon regrowth, the outer tumour radius, R, will increase, while R N will decrease as the necrotic core undergoes decay. Since this parameter regime yields tumour spheroids with proliferating, hypoxic and necrotic compartments prior to radiotherapy, left untreated the tumour will eventually evolve so that the hypoxic annulus will start to re-develop (when c(R N ) = c H ). This sequence of events is summarised in Figure 6. Simulation results reveal the re-emergence of hypoxia following radiotherapy via a fast, transient increase in R H (see Figure 4b). Since hypoxic cells are less radio-sensitive, such a rapid increase in the hypoxic volume could have implications for the response to further doses of radiation. We now characterise this behaviour. The fast initial increase in R H upon regrowth of the tumour spheroid in the Greenspan model after a radiation fraction (see Figure 7) occurs when the necrotic core is undergoing exponential decay and a region of hypoxia is about to develop outside the necrotic core. We analyse this situation by supposing that (without loss of generality at t = 0) c(R N ) = c H . Then solution of Equations (14), (18) and (19) yields while Equation (15) reduces to give and the internal free boundaries R H (t) and R N (t) satisfy with R H (0) = R N (0), and R(0) defined by Equation (26).
With R N (t) defined by Equation (27), we seek approximate solutions for R and R H . Differentiating Equation (26) with respect to t we obtain an ODE for R H (t), which is singular in the limit as R H → R N : We investigate the behaviour in this limit by making the change of variables S = R H R N with S = 1 at t = 0. Then Equations (25) and (28) (30) We construct approximate solutions to Equations (29) and (30) in the limit when S is close to 1. Introducing the small parameter (0 < 1), we consider the short timescale t = 2 τ and propose expansions for S and R in this boundary layer of the form and The choice of timescale can be justified by a dominant balance argument, allowing us to regularise the ODE for S at leading order. Note that dR dτ = o( 2 ), so we deduce that dR 0 dτ = 0 = dR 1 dτ , . We also find that R N has no o( ) term, since expanding Equation (27) gives Turning our attention back to Equation (30), we find at leading order where f (R 0 , R N 0 ) is defined via Equation (35). Integrating Equation (36) gives  (14)- (23) showing the regrowth of the tumour spheroid after radiotherapy and the fast initial increase in R H . We note that during regrowth the mismatch in the oxygen tension on the necrotic core boundary is eventually resolved and the tumour resumes standard Greenspan growth dynamics.
So for t 1, This approximate solution for R H (t) and the numerical solution obtained by solving the full problem are in good agreement (Figure 7). As a result of this analysis we conclude that whenever a hypoxic region re-emerges within a tumour post-radiotherapy, it does so rapidly over a short timescale on which both the outer tumour radius and the radius of the necrotic core do not change dramatically. This phenomenon arises naturally from the model. We note that similar analysis was performed on the original, untreated growth equations by Byrne et al. (Byrne and Chaplain, 1998).

Single hit regrowth
We now consider the effect of irradiating a small tumour spheroid composed entirely of proliferating cells and its subsequent regrowth. This situation is relevant for tumours of radius R such that 0 < R 2 < 6 Γ (c ∞ − c H ) and R H = R N = 0, since when R = 6 Γ (c ∞ − c H ), c(0) = c H and therefore larger tumours will contain a hypoxic region and/or necrotic core. Solving Equation (14) subject to boundary conditions (18)-(19) we deduce that the oxygen profile for this tumour composition is given by Substituting from Equation (39) into Equation (15), with R H = R N = 0, we arrive at the following ODE for R(t): with solution 15 . We can use Equation (41) to calculate the time taken for a tumour spheroid of radius 0 < R 2 < 6 Γ (c ∞ − c H ) following a single fraction of radiation of dose d to regrow to its initial (i.e. pre-radiotherapy) size. The tumour radius immediately after irradiation is given by γ(d)R, where γ(d) = exp(− 1 3 (αd + βd 2 )) is the survival fraction given by the linear quadratic formula. It is straightforward to show that the time, ∆t, until the tumour regrows to its original size is given by where γ = γ(d). We assume that the parameters are such that the tumour spheroid was growing pre-irradiation and therefore R 2 < 15 Γ (c ∞ − λ A ) (from (40)). In this parameter regime, the logarithm in Equation (42) is defined and ∆t > 0. By extension, we require that c ∞ > λ A since otherwise the tumour spheroid shrinks for all values ofR and a viable tumour of any size cannot be supported. We note for future reference that Differentiating Equation (42) with respect toR, we obtain Therefore, at least for smallR, d∆t dR γ(d) > 0, and so the regrowth time, ∆t, is an increasing function ofR for fixed dose d. This holds for all 0 <R 2 < 6 Γ (c ∞ − c H ) if the inequality 3c ∞ − 5λ A +2c H > 0 is satisfied. Then Equation (42) represents an increasing family of curves where the minimum bounding curve is given by Equation (43). That is, for 0 <R 1 <R 2 < 6 Γ (c ∞ − c H ) and fixed dose d, ∆t(d; R 1 ) < ∆t(d; R 2 ).
Returning to Equation (42), we see that for an initial radiusR, a strategy that combines a dose d with an inter-fraction time less than ∆t will result initially in a net shrinkage of the tumour throughout treatment (i.e. the tumour volume is smaller at the delivery of each radiation fraction). Conversely, if we wait longer than ∆t to re-irradiate then the tumour will have grown larger than its original size. Thus, Equation (42) defines a curve in the (d, ∆t)-plane for treatment protocols which give rise to periodic behaviour for a given initial tumour radiusR, since the cell kill induced by radiotherapy is exactly balanced by the regrowth between fractions. This curve is indicated in Figure 8. We note that the curve describing protocols that yield periodic behaviour becomes concave for large R, and in particular, as R tends to its steady state value, R * say, then for a given dose, d > 0, ∆t(d) → ∞.

Periodic surface in 'treatment space'
In the previous section we investigated the response of tumours of fixed initial radiusR to different fractionation protocols as specified by a dose, d, and inter-fraction time, ∆t. The (d, ∆t)-plane can be partitioned into regions in which the tumour volume increases or decreases during treatment. However, we note that the location of these regions depends on the tumour radius at the time of irradiation,R. As such, when considering the behaviour of a tumour throughout an entire course of radiotherapy, it is not a line in the (d, ∆t)-plane that we are interested in, but a surface in 'treatment space', T ⊂ R 3 , with components (R, d, ∆t). For a protocol delivering n fractions of radiation, we can identify the response of a given tumour with a discrete trajectory, (R(t i ), d i , ∆t i ) ∈ T, through treatment space. Note, we consider the trajectory as discrete points pre-irradiation so that R(t i ) = R(t i ) − in our previous notation. Here t i , d i and ∆t i are the time, dose and inter-fraction time, respectively, of the ith fraction, for i = 1, ..., n. The union of the curves in the (d, ∆t)-plane determining periodicity define a surface, S, in T (see Figure 9). Now, as a tumour progresses through the course of radiotherapy, points on its trajectory that lie above this surface indicate net growth of the tumour by the time of delivery of the next radiation dose, whilst points below the surface result in a more desirable decrease in tumour volume. Therefore, for a given tumour spheroid, S defines a surface that partitions T into treatment protocols that cannot halt tumour progression and those that lead to tumor decay.
With this understanding more general statements about the outcome of different fractionation protocols can be made. We first consider dosing schedules in which the same dose d is administered at constant intervals ∆t. In this case, if a point in T lies above S, then the tumour spheroid will grow until it converges to a point on S, at which time it becomes periodic ( Figure 9). Below S, then the outcome of radiotherapy depends on the boundary curve of S as R → 0 (see Equation (43)). If the fractionation schedule lies above this curve in the (d, ∆t)-plane then the system will eventually converge to the corresponding point on S. However, d and ∆t below this can be chosen such that the tumour decays to arbitrarily small volumes (Figure 9). We note that Equation (43) still holds for this analysis even though its derivation requires an initial tumour composed entirely of viable cells. This is the case since, in the model, necrotic material is only created via nutrient starvation and not as a result of irradiation. As such, for ongoing treatment protocols that result in a decreasing tumour volume, as t → ∞, R R N . We can extend this concept to more general treatment schedules. If there exists m ∈ N such that for all i > m the point (d i , ∆t i ) of the ith fraction lies below the curve given by Equation (43), then treatment is sufficient to drive the tumour volume arbitrarily close to 0. However, using Figure 4b as an illustrative example, we observe a simulation in which the weekend break of the standard fractionation protocol corresponding to ∆t = 3 days gives rise to the periodic orbit.
When predicting how a tumour may respond to radiotherapy, two quantities of interest are the potential volume doubling time, T pot , and the survival fraction after 2 Gy of radiation, SF 2 , of the tumour. We translate these definitions into the model and identify tumour characteristics Figure 9: Simulated trajectories (green lines) in T for two different constant dosing schedules (fractions plotted as red points) applied to the same tumour spheroid (i.e. same model parameters). A dose of 3 Gy every 3 days (top trajectory) results in a tumour that grows until it reaches a periodic state, while a dose of 2 Gy delivered daily (bottom trajectory) is sufficient to drive the tumour volume arbitrarily close to zero. Surface describes treatment protocols which give rise to periodic behaviour. for which this model would eventually predict regrowth over the weekend.
The survival fraction at 2 Gy, SF 2 , for normoxic tumour cells is given by the linear-quadratic formulation. In order to find an expression for T pot , we consider the tumour dynamics for small tumour volumes that are in an 'exponential growth' phase. For small R, which yields a volume doubling time of Substituting these expressions into Equation (43) for a standard 2 Gy dose we obtain We note that this expression holds for tumour growth models in which the growth of small tumours is exponential and radiation-induced cell death is modelled as an instantaneous volume loss.
If we now consider Equation (45) in the context of the weekend break (∆t = 3) within the standard fractionation protocol, then we see that −T pot log(S 2Gy ) < 3log2 defines a class of tumour characteristics for which regrowth of the tumour over the weekends impedes tumor decay ( Figure  10).

Discussion
In vivo tumours are highly heterogeneous cellular entities characterised by high inter-patient variability. Allied to this, the local efficacy of radiotherapy delivered to the tumour site is affected by a number of variables associated with the tumour's heterogeneous composition and microenvironment. In particular, the local oxygen concentration can significantly influence radiationinduced cell death, with well-oxygenated regions being shown to exhibit up to threefold greater Figure 10: Region for which tumours will eventually exhibit net growth over the weekend (shown in red). radiosensitivity than hypoxic tumour populations. With this in mind, in this paper we have presented a simple, spatially-resolved model in order to investigate the effects of tumour composition on radiotherapy response. The discussed model builds on the tumour growth model first proposed by Greenspan (Greenspan, 1972). We extend this to incorporate the effects of radiotherapy taking into account spatially-varying radiosensitivities.
Numerical simulations and mathematical analysis of the model reveal how the tumour's growth dynamics and spatial composition change throughout treatment. Heterogeneity within the tumour not only affects the initial response to radiotherapy, but also how this response changes throughout the duration of the treatment protocol (see Figure 5b). For parameter regimes in which hypoxia reemerges during treatment, a rapid transient increase in the width of the hypoxic annulus that is observed. This behaviour arises naturally from the model and the underlying process driving this phenomenon remains to be elucidated. The model more generally also classifies protocols that may result in tumour progression, a non-zero periodic tumour volume, or tumor decay. We identify a surface in 'treatment space' dependent on tumour-specific growth and radiosensitivity parameters and determine that successful protocols correspond to those that remain below this surface throughout treatment. The wide variety of dynamics observed suggests that spatial heterogeneity may be important for simulating tumour response to radiotherapy and, in particular, for making clinical predictions.
The model presented in this paper makes numerous assumptions and simplifications about the underlying biology. For radiation-induced cell death we take the common approach of modelling this process as an instantaneous effect. However, the linear-quadratic model was established to determine long-term clonogenic survival after radiotherapy. Biologically, cell death after radiation may occur via a number of different mechanisms, with many irradiated cells dying only after attempting mitosis one or more times (Joiner and van der Kogel (2009), Chapter 3). Consequently, radiation-induced cell death may not elicit the instantaneous volume reduction modelled here. In future work we will model this process in more detail in order to describe the short term response to radiotherapy and the corresponding spatial changes in tumour composition more accurately. Additionally, a more detailed model will allow us to relax the restrictions of the spherically symmetric geometry associated with the current model, making it more applicable to modelling vascular tumour growth and in vivo responses to radiotherapy.
However, for a model to be useful in making patient-specific predictions in the clinic, parameters must be identifiable with respect to the limited clinical data typically available. As such, current research in this area often focusses on phenomenological ODE models with few parameters to be estimated. Typically these models contain no information about the spatial heterogeneity in tumour composition and radiotherapy response. With this in mind, we also propose a compar-ison of the developed spatially-resolved model with these phenomenological approaches as future work. We aim to identify situations in which the spatially-resolved and spatially-averaged models agree well, and those in which there is a significant difference. The tumour composition changes observed in model simulations suggests that averaged parameter values in simple, phenomenological models may not sufficiently capture the tumour dynamics during treatment for some tumour compositions and parameter combinations.
While more sophisticated models may be difficult to parametrise in practice, they have the potential to increase biological insight and inform further modelling studies. Spatially-resolved models, such as the one presented in this paper and the future work proposed to generalise some of the simplifying assumptions made here, may aid in the development of alternative clinicallyfocussed models which capture more of the key features than existing phenomenological models. More complex models incorporating more biological detail may also be used to generate data for in silico testing of ODE models and model selection paradigms; comparing the quality of fit and future predictions of a range of simple models against a known 'ground-truth'.

Acknowledgements
This work was supported by the Engineering and Physical Sciences Research Council [grant number EP/G037280/1]. TL would also like to thank the Moffitt Cancer Center, where some of this work was undertaken, for their hospitality. rate, Γ, and that cells in the necrotic core do not consume oxygen. Since oxygen diffuses on a much shorter timescale than tumour growth we assume that the oxygen concentration is in a quasi-steady state. Hence, c(r, t) satisfies (46), where H(.) is the Heaviside function, with associated boundary conditions (50) and (51) imposing symmetry at r = 0 and a Dirichlet boundary condition at r = R(t).
We assume that rates of cell proliferation and death within the tumour are determined by the local oxygen concentration. Proliferation occurs at a rate proportional to c where there is sufficient oxygen supply (c > c H ). Both apoptosis and necrosis contribute to cell death within the tumour, with apoptosis occurring at a constant rate throughout the tumour and necrosis localised to the necrotic core (where c < c N ). Degradation of the necrotic core is assumed to result in material that is freely permeable throughout the tumour spheroid. We assume that adhesion and surface tension forces acting on the tumour cells maintain the shape of the tumour spheroid, and that these same forces push cells inwards to compensate for the outward flux of necrotic material from the necrotic core. The evolution of R(t) is then given by the mass balance in Equation (47), with λ A , λ N > 0, accompanied by the initial condition (52).
The radii at which the tumour becomes hypoxic and necrotic, R H and R N , respectively, are determined as contours of the oxygen concentration profile (Equations (48) and (49)). For situations in which the tumour spheroid is small enough that one or both of these contours do not exist, we define the corresponding radius to be 0.
An analytic solution can be found for c(r, t) that depends on R, R H and R N , and so this system can be reduced to an ODE for R and algebraic equations for R H and R N . The resulting equations in the case of a fully developed, 3-layer tumour spheroid are given in Equations (53)-(57). The system of equations for the earlier growth phases are similar and can be obtained in the same manner.  Table 3. In using partial pressures of oxygen rather than concentration, we follow Grimes et al. (Grimes et al, 2014) and use Henry's law to convert between the two, so that p = Ωc, with Ω = 3.0318x10 7 mmHgkgm −3 . The oxygen concentration profile at each time point is a monotonic function increasing outwards from the tumour centre. Initially, for very small, avascular tumour spheroids, the entire tumour is made up of viable, proliferating cells (case A). We notice that as the tumour grows and so the oxygen concentration at the centre of the tumour decreases, a hypoxic region within the tumour begins to develop (case B) followed by a central necrotic core when the oxygen concentration falls so low as to be unable to support viable cells (case C). These different regions, and their varying responses to radiation, require careful consideration when extending Greenspan's model to account for radiotherapy. Further discussion of the different phases of growth is given in (Byrne, 2012). Table 3: Parameter values for Greenspan's model of tumour spheroid growth. (Note: dimensional apoptosis and necrosis rates, s −1 , are given by sλ A and sλ N , respectively. p ∞ , p H and p N are partial pressures which, for consistency, are then converted into concentrations c ∞ , c H and c N , respectively. For further details see Appendix A.1.)

A.2 Parameter sweep of Greenspan growth dynamics
We initially investigate the sensitivity of the standard Greenspan growth dynamics to some of its key parameters. In particular, as we sweep over a region of parameter space, we observe the resulting long-time, steady-state properties of the tumour spheroids. Specifically, we investigate how the availability of oxygen, c ∞ , the oxygen consumption rate of the cells, Γ, and the rates of apoptosis and necrosis, λ A and λ N , respectively, affect the growth of the tumour under this model. We choose appropriate intervals for each parameter, in each case encompassing the corresponding value shown in Table 3, and systematically explore the resulting region of parameter space, solving the model numerically. We plot heatmaps in parameter space in order to observe the behaviour of a variety of quantities of interest of the tumours at steady state (Figures 11 & 12). We have omitted here the results of varying both c ∞ and λ N since, in the case of c ∞ , increasing the parameter value simply has the effect of an increase in overall tumour size and the width of the corresponding proliferating rim as might be expected, while, for the range of values swept over, λ N had little effect on the observed tumour characteristics.
Broadly speaking, we see in Figure 11 that relatively increasing the amount of oxygen available to the tumour cells (by decreasing the consumption rate Γ), or decreasing the death rate, λ A , results in larger steady state tumours, as we might expect. Point A marks the location of the parameter values from Table 3. Extreme points in the region of concern in parameter space are labelled B, C and D and the corresponding tumour evolution profiles shown. We see that the greater relative availability of oxygen and the low death rate at point B results in a large, fullydeveloped tumour spheroid, whereas the high rate of apoptosis at D is such that a viable, tumour cell population cannot be sustained.
For a given Γ, we can determine the onset of necrosis as the radiusR = 6 Γ (c ∞ − c N ). At the point C, and any other point on a vertical line through C, the rate of apoptosis is such that the tumour reaches a steady state of radiusR, and is given by λ A = 1 5 (3c ∞ +2c N ). Beyond C, at lower death rates (eg B), we enter a new phase of growth that includes hypoxia and necrosis. The rapid, transient increase in R H /R N at onset of hypoxia/necrosis (c.f. asymptotics in Section 3.2) means that for each tumour the width of the proliferating rim decreases until it reaches a steady state and as such C represents a 'local maximum'. Figure 12 shows the corresponding heatmaps for both the steady state volumes and volume fractions for each constituent of the tumour spheroid.