Investigating the Influence of Growth Arrest Mechanisms on Tumour Responses to Radiotherapy

Cancer is a heterogeneous disease and tumours of the same type can differ greatly at the genetic and phenotypic levels. Understanding how these differences impact sensitivity to treatment is an essential step towards patient-specific treatment design. In this paper, we investigate how two different mechanisms for growth control may affect tumour cell responses to fractionated radiotherapy (RT) by extending an existing ordinary differential equation model of tumour growth. In the absence of treatment, this model distinguishes between growth arrest due to nutrient insufficiency and competition for space and exhibits three growth regimes: nutrient limited, space limited (SL) and bistable (BS), where both mechanisms for growth arrest coexist. We study the effect of RT for tumours in each regime, finding that tumours in the SL regime typically respond best to RT, while tumours in the BS regime typically respond worst to RT. For tumours in each regime, we also identify the biological processes that may explain positive and negative treatment outcomes and the dosing regimen which maximises the reduction in tumour burden.


Introduction
Understanding the biological mechanisms underpinning cancer and developing effective therapeutic protocols to improve patient prognosis are fundamental aims of cancer research. Existing treatment modalities, such as radiotherapy (RT) and chemotherapy (CT), are applied via highly-regulated dosing protocols (National Cancer Institute 2019 to avoid damaging healthy tissue, while maximising treatment effect. Nonetheless, the efficacy of both RT and CT is limited by their intolerable side-effects. Further, inter-tumour heterogeneity can significantly influence sensitivity to treatment. Investigating how different growth mechanisms may affect response to treatment is, therefore, an important step towards overcoming barriers to treatment efficacy. In this paper, we investigate how two distinct growth-rate limiting mechanisms, namely growth arrest in response to nutrient insufficiency and to competition for space, impact tumour response to RT. A dynamic model of tumour growth that distinguishes between mechanisms of tumour control. Regardless of their biological complexity, existing models of solid tumour growth typically describe a single mechanism by which a tumour may reach an equilibrium. For example, the models developed by Enderling et al. (2006), Hahnfeldt et al. (1999), Liu et al. (2021), Milzman et al. (2021) and Zahid et al. (2021) predict growth arrest due to a cessation of proliferation (with no explicit cell death), while those proposed by Drasdo and Höhme (2005), Greenspan (1972), Hillen et al. (2013) and Lewin et al. (2020) predict growth arrest due to the balance of cell proliferation and cell death.
In previous work (Colson et al. 2022), we developed a model of solid tumour growth that simultaneously describes growth arrest due to nutrient inhibition, when the net rates of cell proliferation and death are equal (and strictly positive), and growth arrest due to contact inhibition, when the net rate of cell proliferation becomes zero, with no cell death. We assumed that the system is well-mixed and neglected angiogenesis and vascular remodelling. In particular, we viewed the vascular volume as a parameter which influences nutrient and space availability and, therefore, the tumour's carrying capacity. As such, the model does not capture the co-evolution of the vasculature with tumour cells that is observed in vivo and that can allow nutrient limited tumours to continue to grow (see Discussion). Nor does it describe how secondary tumours are established, i.e., metastasis.
The model comprises two time-dependent ordinary differential equations (ODEs) for the tumour volume, T (t), and the oxygen concentration, c(t), and can be written as follows: rate of tumour cell proliferation where Denoting the total available space by S max (m 3 ) and the vascular volume by V 0 (m 3 ), the rate of tumour cell proliferation is assumed to be proportional to the available free space, S max − T − V 0 , and to the oxygen concentration, c, with proportionality constant q * 2 (kg −1 min −1 ). If c drops below a threshold value, c * min (kg m −3 ), then cells die at a rate proportional to the difference between c and c * min , with proportionality constant δ * 1 (m 3 kg −1 min −1 ). Further, oxygen is supplied to the tumour at a rate proportional to V 0 and the difference between the oxygen concentration in the vasculature, c * max (kg m −3 ), and in the tumour. The parameter g * (m −3 min −1 ) is the rate of oxygen exchange per unit volume area of blood vessel. Finally, oxygen is consumed by tumour cells for maintenance at a rate proportional to c, with rate constant q * 1 (m −3 min −1 ), and for proliferation at a rate proportional to the proliferation rate, with conversion factor k > 0 defined such that q * 3 = q * 2 /k (m −6 min −1 ). Since the model (1), (2) distinguishes between two mechanisms for growth-control, it can be used to investigate how they impact tumour response to treatment. Therefore, in this work, we extend Eqs. (1), (2) to account for the biological effects of RT. Radiobiology. RT is used to treat more than 50% of cancer patients (Maier et al. 2016). It involves the delivery of energy rays, via small doses called fractions, at regular time intervals and over a fixed period of time, to a region of the body comprising both cancerous and healthy tissue. Radiation protocols are, therefore, designed to balance treatment efficacy and undesirable side-effects in normal tissues. While a conventional fractionation schedule consists of 2 Gy doses delivered Monday to Friday for up to 7 weeks (Ahmed et al. 2014), the dose, dosing frequency and treatment duration can be varied to deliver a fixed total dose. The latter is often termed the Maximum Tolerated Dose, i.e., the highest dose which does not cause adverse side-effects (Gad 2014).
Radiation induces direct and indirect cytotoxic effects by causing DNA damage to cancer cells that is fatal if left unrepaired. Direct effects arise from interactions between ionising particles and DNA and indirect ones from interactions between ionising particles and water, which create reactive oxygen species that subsequently react with DNA. Indirect effects are the most common, which is why hypoxic, i.e., poorly oxygenated, tumours are often radio-resistant (Graham and Unger 2018). Intratumoural oxygen levels are a key factor influencing tumour radio-sensitivity, and another key factor are tumour cell proliferation rates, as cells that are in the G2 or mitosis phases of the cell cycle are the most sensitive to RT.
RT can also affect the tumour vasculature, with increases in angiogenesis observed at low doses (Marques et al. 2020) and vascular damage and necrosis observed at high doses (Stolz et al. 2022;Venkatesulu et al. 2018). In this work, we neglect the effect of RT on vasculature in order to focus on evaluating how nutrient-and contact-inhibited growth control affect the sensitivity of tumour cells to treatment with RT. Mathematical modelling of tumour response to radiotherapy. A number of mathematical models have been proposed to describe tumour response to RT. Key aims of these modelling efforts include studying specific RT protocols (Enderling et al. 2010;Jeong et al. 2017;Lewin et al. 2018;Prokopiou et al. 2015;Rockne et al. 2009), designing patient-specific RT dosing schedules (Alfonso et al. 2014;Belfatto et al. 2018) and investigating the influence of inter-and/or intra-tumour heterogeneity on tumour sensitivity to RT (Alfonso and Berk 2019;Celora et al. 2023;Powathil et al. 2012;Enderling et al. 2009;Watanabe et al. 2016).
While the purpose of these modelling approaches may differ, they are all based on the common assumption that RT inflicts instantaneous cell death on tumour cells and the cell kill is modelled using the Linear-Quadratic (LQ) model (McMahon 2018). The LQ model states that the fraction, S L Q , of (tumour) cells that survive exposure to a dose D (Gy) of radiation is given by where α ≥ 0 and β ≥ 0 are tissue-specific radio-sensitivity parameters. These parameters are typically derived from cell survival data collected at a small number of time points in in vitro two-dimensional (2D) monolayer or three-dimensional (3D) spheroid experiments. As such, they provide information about the long-term proportion of cell death rather than how the cell death rate changes over time. In contrast, a time-dependent description of RT cell kill can account for different types of damage (direct vs. indirect), damage repair and cell death following insufficient repair (Curtis 1986;Goodhead 1985;Neira et al. 2020;Tobias 1985). Such a description facilitates the study of the evolution of tumour composition during treatment, as we may keep track of changes in healthy, damaged and dead tumour cell populations. In this paper, we follow Neira et al. (2020) and adopt a time-dependent description of RT. Paper structure. This paper is structured as follows. In Sect. 2, we extend the tumour growth model defined by Eqs. (1), (2) to account for the biological effects of RT. We summarise the key features of the model dynamics in the absence of treatment in Sect. 3. Then, we investigate the response of tumours characterised by different growthlimiting mechanisms in Sect. 4, initially performing a numerical study of tumour response during RT and, subsequently, looking at post-treatment growth dynamics via a steady state analysis and complementary numerical study. The paper concludes in Sect. 5, where we discuss our findings and outline possible avenues for future work.

The Mathematical Model
In this section, we incorporate tumour response to RT into the growth model (1), (2). We follow the approach outlined in Neira et al. (2020) and adopt a time-dependent description of radiotherapy. In more detail, we introduce the dependent variables T S and T R to denote tumour cells that have been, respectively, sub-lethally and lethally damaged by RT. We suppose that the tumour is exposed to a total dose D (Gy) of RT at a constant rate R over the time period t R ≤ t ≤ t R + δ R (min) so that Let = T + T S + T R + V 0 be the total tumour volume. We propose the following system of time-dependent ODEs to describe tumour growth and response to RT (see also the schematic in Fig. 1): where H is the Heaviside function defined in (3). We assume that undamaged tumour cells, T , proliferate, die due to nutrient insufficiency and consume oxygen for proliferation and maintenance as in Eqs. (1), (2). They suffer sub-lethal and lethal damage during irradiation at rates proportional to the oxygen concentration, c, and the RT dose rate, R, with proportionality constants ν * > 0 and λ * > 0, respectively. We further suppose that sub-lethal damage is either repaired at a constant rate μ * > 0, or leads to tumour cell death via two distinct pathways. First, sub-lethal damage may become lethal as it accumulates at a rate proportional to the oxygen concentration, c, and the RT dose rate, R, with proportionality constant λ * S > 0. Second, sub-lethally damaged cells, T S , may also undergo mitotic catastrophe (MC) if they attempt to divide with mis-or un-repaired DNA damage; we assume this occurs at a constant rate ξ * > 0.
Sub-lethally damaged cells, T S , also consume oxygen for maintenance and proliferation, proliferate and die due to nutrient insufficiency similarly to undamaged cells, T , although at different rates. More specifically, they proliferate at a rate proportional to the oxygen concentration, c, and the available space, with proportionality constant q * 2,S = θ 2 q * 2 > 0, with θ 2 ∈ (0, 1). The latter ensures that damaged cells proliferate more slowly than undamaged cells as they expend more energy repairing RT damage Fig. 1 Schematic showing the interactions between undamaged, damaged and dead tumour cells, T , T S , T R , respectively, in response to RT and the proliferation of T and T S cells as described in the model defined by Eqs. (5)-(9). R denotes the RT dose rate defined by (5) and c denotes the intratumoural oxygen concentration (Colour figure online) than proliferating. Accordingly, sub-lethally damaged cells consume oxygen for maintenance at a rate proportional to c, with rate constant q * 1,S = θ 1 q * 1 > 0, where θ 1 > 1 as these cells require more energy to repair RT damage. They also consume oxygen for proliferation at a rate proportional to the proliferation rate, with conversion factor Here, we assume the same conversion factor for T and T S cells, for simplicity. Since q * 2,S = θ 2 q * 2 and q * 3 = q * 2 k , we also have q * 3,S = θ 2 q * 3 , i.e., damaged cells consume less oxygen for proliferation than undamaged cells. Lastly, as for T cells, T S cells die from nutrient insufficiency when c < c * min , at a rate proportional to the difference between c and c * min , with proportionality constant δ * 1,S > 0. Lethally-damaged cells, T R , are considered to be dead: their damage cannot be repaired, they do not consume oxygen or proliferate, but they occupy space and are degraded at a constant rate η * R > 0. One final and important assumption we make is that radiation only affects tumour cells, i.e., we neglect any effects RT may have on tumour angiogenesis, vascular remodelling and injury. This simplifying assumption enables us to focus on elucidating how mechanisms of growth arrest influence one particular type of tumour response to RT, i.e., the cellular response.

Non-dimensionalisation
We non-dimensionalise Eqs. (5)-(9) by introducing the following scalings: The timescale of interest is the timescale for the duration of RT (τ = 1 min). We choose this timescale since we seek to describe the damage and repair associated with RT. The maximum dose rate is also fixed to be R max = 1 Gy/min (Konopacka et al. 2016). Then, given that q * 1,S = θ 1 q * 1 , q * 2,S = θ 2 q * 2 and q * 3,S = θ 2 q * 3 and dropping hats for notational convenience, we obtain the following dimensionless system: where we have introduced the following dimensionless parameter groupings:

Defining the Dimensionless Model Parameters
This paper focusses on studying the impact of two distinct growth arrest mechanisms on the qualitative tumour response to RT. We, therefore, fix parameters related to tumour cell responses to RT at the values stated in Table 1. The values of ν, λ, λ S , μ and η R are based on values found in the literature (Neira et al. 2020;Steel et al. 1987) and we assume λ = λ S , for simplicity. We also set ξ = 5 × 10 −4 so that cells that undergo mitotic catastrophe have a half-life of approximately 24h. Here, we implicitly assume that the average duration of the cell cycle in healthy cells (Bernard and Herzel 2006) and cancer cells are approximately the same. The parameters that define the RT dosing schedules (e.g., the dose rate, R) are summarised in Sect. 4.1. Further, we define tumour growth parameters as in our previous work (Colson et al. 2022). In particular, we fix the values of c min , g and k based on parameter values found in the literature and preliminary numerical simulations. For q 1 , q 3 and V 0 , we consider a biologically realistic range of possible values, motivated by arguments outlined in Colson et al. (2022), as these parameters have a strong influence on the qualitative behaviour of the model in the absence of treatment (see Sect. 3). We make the additional simplifying assumption that, as for the undamaged cells, δ 1,S = q 2,S . Finally, we set θ 1 = 10 and θ 2 = 0.1 to represent a 10-fold increase in oxygen consumption for maintenance and a 10-fold decrease in oxygen consumption for proliferation in damaged cells.  Fig. 2 For q 1 ∈ {0.1, 0.5, 1}, we show the regions of (V 0 , q 3 )-space in which only the stable NL steady state exists (blue), only the stable SL steady state exists (red) and both the stable NL and SL steady states co-exist (purple). The solid and dashed red lines represent the boundaries between the three regions. For Tumours defined by parameter sets A 1 -C 1 have values of q 3 which are sufficiently larger than q 1 so that there is bi-stability, while, for A 2 -C 2 , bistability does not occur (Colour figure online) et al. (2022), we showed that this model admits two non-trivial, stable steady states (see Appendix A): (i) a nutrient limited (NL) steady state, attained when cell proliferation balances cell death due to nutrient starvation; (ii) a space limited (SL) steady state, attained when cell proliferation ceases due to lack of space, with no cell death. For these solutions to be physically realistic and lie in the appropriate nutrient regime, they must satisfy 0 ≤ T ≤ 1 − V 0 and either 0 ≤ c < c min (for the NL steady state) or c ≥ c min (for the SL steady state). Imposing these conditions, we find that admissible NL and SL steady state solutions lie in different regions of parameter space. These regions are defined by the values of q 1 , q 3 and V 0 , with no qualitative changes when c min and g vary (Colson et al. 2022). Figure 2 depicts these regions in (q 3 , V 0 )-space for three values of q 1 and fixed values of c min = 0.01 and g = 5.
Given q 1 , there exists a threshold value of V 0 , which we denote by V N , which is independent of q 3 , such that only NL steady states exist for 0 < V 0 ≤ V N . Tumours in this region of parameter space (e.g., A 1 and A 2 ) are said to be in a NL growth regime. Further, given q 1 , and for q 3 sufficiently large relative to q 1 , there exists another threshold value of V 0 , which we denote by V S , such that only SL steady states exist for V 0 ≥ V S > V N . Tumours in this region (e.g., C 1 , B 2 and C 2 ) are said to be in a SL growth regime. Further, where V N < V 0 < V S , the NL and SL steady states co-exist. A tumour in this region (e.g., B 1 ) may evolve to either steady state, depending on its initial conditions. We consider such tumours to be in a bistable (BS) growth regime. Finally, for q 3 q 1 , we have that, for V 0 > V N , a unique steady state exists and it is of SL type. Figure 3 shows the time evolution of the tumour cell volume, T , and the intratumoural oxygen concentration, c, and the corresponding bifurcation diagrams for tumours A 1 -C 1 and A 2 -C 2 . In all cases, the NL steady state values for T and c (T * and c * , respectively) are smaller than the SL ones. This is consistent with the assumption that, in the absence of angiogenesis, well-oxygenated tumours attain larger volumes than poorly-oxygenated tumours. Further, tumours in a BS regime evolve to their NL steady state for initial conditions satisfying 0 < T (0) T * , which we use to simulate tumour growth. As the values of T * and c * for NL tumours increase with V 0 , tumours (0)) = (0.05, 0, 0, 1) and plot the evolution of the tumour volume and oxygen concentration in time. In (a), (q 1 , q 3 , V 0 ) correspond to points A 1 -C 1 in Fig. 2 and, in (b), they correspond to points A 2 -C 2 . The remaining model parameters are fixed at the default values in Table 1.The bifurcation diagrams show how the steady state values for T and c change as V 0 varies for (q 1 , q 3 ) corresponding to A 1 -C 1 in (a) and for (q 1 , q 3 ) corresponding to A 2 -C 2 in (b). In both cases, the NL steady state increases with V 0 and is smaller than the SL steady state, which decreases with V 0 . The tumour in a BS regime (B 1 ) grows to its NL steady state (Colour figure online) in a BS regime will grow to larger volumes than tumours in a NL regime (and smaller volumes than tumours in a SL regime). We also note that, in BS regimes, there is a large jump in T * and c * at V S , the threshold value of V 0 separating the BS and SL regimes. In contrast, in monostable regimes, T * and c * depend continuously on V 0 .
In summary, in the absence of treatment (i.e., R ≡ 0), Eqs. (10)-(13) describe three possible scenarios for tumour growth: (i) nutrient limited growth, where the tumour grows to a NL steady state at which proliferation balances cell death due to nutrient insufficiency; (ii) space limited growth, where the tumour grows to a SL steady state at which proliferation ceases due to space constraints; (iii) bistable growth, where a tumour grows to a NL steady state given physically realistic initial conditions (0 < T (0) T * ). In vivo, we expect that tumours growing in poorly-perfused regions (e.g., breast, bone) and/or that elicit a weaker angiogenic response would be in the NL regime, while tumours growing in well-perfused regions (e.g., lung, liver, brain) and/or that elicit a stronger angiogenic response would be in the SL regime. Tumours growing in well-perfused regions that are highly proliferative and, therefore, more likely to outgrow their nutrient supply, would be in the BS regime. In Sect. 4, we investigate how tumours in these three growth regimes respond to RT.

Methods
Our aim is to understand the qualitative response of tumours in nutrient limited (NL), space limited (SL) and bistable (BS) regimes to a range of fractionated radiotherapy (RT) treatments. As a first step, we create three virtual tumour populations as follows. We first fix all tumour growth model parameters, except q 1 , q 3 and V 0 , at the default values stated in Table 1. We then also fix V 0 = 0.0005, V 0 = 0.005 and V 0 = 0.00275 for tumours in the virtual NL, SL and BS regimes, respectively. Allowing q 1 and q 3 to vary, we generate three virtual tumour populations of size N = 250 by randomly selecting N (q 1 , q 3 ) pairs, corresponding to the NL, SL and BS regimes, respectively. Each pair is formed by sampling pairs of values of q 1 and q 3 from uniform distributions whose lower and upper bounds depend on the growth regime (NL, SL or BS), the value of V 0 and, for sampling q 3 in the BS and SL regimes, the value of q 1 (see Appendix B).
We then define the RT protocols of interest. We vary the dose amount D ∈ 0, 5 (Gy) and the number of doses per week N f rac = {1, 3, 5}. We assume that each dose is administered in δ R = 10 min and, therefore, we vary the dimensionless dose rate, We also suppose that all fractions are applied at the same time of day and the first weekly fraction is applied on Mondays, with subsequent fractions applied at equally spaced time intervals during Monday to Friday (e.g., 3 doses per week corresponds to doses on Monday, Wednesday and Friday). Further, the duration of each fractionation schedule is determined so that the total dose administered is 80 Gy (or the closest multiple of D to 80 Gy).
For each set of tumours and each RT protocol, we solve Eqs. (10)-(13) numerically for t ∈ (0, t * ], t * > 0, using ODE45, a single step MATLAB built-in solver for non-stiff ODEs that is based on an explicit Runge-Kutta (4,5) formula, the Dormand-Prince pair (Dormand and Prince 1980). (The code is available at: https://github. com/chloeacolson/InvestigatingTumourResponsesRadiotherapy.) For simplicity, we impose the initial conditions where T * and c * are, respectively, the steady state tumour volume and oxygen concentration in the absence of treatment. All RT parameters are fixed at the default values listed in Table 1.
For each simulation, we recordT ,T S andT R , the mean undamaged, damaged and dead cell volumes in the last week of treatment, respectively. We also define the percent change in (mean) viable and total cell volumes between the start and the end of treatment as follows We also quantify the end-of-treatment tumour composition (relative to the total tumour volume at the start of treatment) as follows We note that the variables defined in (17) can be used to describe total = (%T + %T S + %T R + %V 0 ) − 100. Henceforth, for brevity, we refer to the relative changes in tumour cell volumes, as defined in (16) and (17), as changes in tumour cell volumes. Finally, we recordc, the mean oxygen concentration in the last week of treatment, and the post-treatment steady state values of all the dependent variables.

Characterising Tumour Response to Fractionated RT
In this section, we investigate the response of tumours in the NL, SL and BS virtual populations to fractionated RT. For each regime, we initially study tumour response to a conventional fractionation schedule consisting of 5 × 2 Gy fractions per week for 8 weeks. In particular, we determine the average response and explore how certain values of q 1 , q 3 and V 0 generate extremal behaviour. We also study the impact of the dose and dosing frequency on tumour response. We consider monostable regimes before looking at the bistable regime.

Tumours in Monostable Regimes: The NL and SL Virtual Tumour Populations
Typical responses to a conventional fractionation schedule. Fig. 4 shows the response of two NL and SL tumours to RT, both of which experience a decrease in viable tumour cell volume, T + T S , during treatment. Since, in both cases, the dependent variables evolve to time periodic solutions within 5 weeks of treatment, we deduce that there is a maximal reduction in the viable cell volume that can be achieved with this fractionation schedule. This maximum reduction, which we quantify using viable , is significantly larger for the SL tumour at approximately 37.6% than for the NL tumour at approximately 4.36%. RT is more effective for the SL tumour as it is better oxygenated than the NL tumour and, hence, it experiences a higher rate of RT cell kill and greater accumulation of dead material, T R , than the NL tumour. Figure 4 also shows that, for both tumours, the oxygen concentration and the viable tumour cell volume decrease when RT is applied. This is because T and T S cells consume oxygen at different rates: we recall that the oxygen consumption rates of sub-lethally damaged cells satisfy q 1,S = 10q 1 and q 3,S = 0.1q 3 . Therefore, changes in tumour composition during treatment will alter the overall oxygen consumption rate of viable tumour cells, leading to transient, or persistent, increases or decreases in the oxygen concentration depending on the values of q 1 and q 3 . Figure 5 shows the distributions of viable and total , following a conventional fractionation schedule, across the NL and SL virtual populations. We note that the behaviour shown in Fig. 4 for specific NL and SL tumours is representative of the average behaviour of each virtual population. In particular, tumours in the SL cohort typically respond well to treatment, with median and (Q1, Q3) values of viable equal to −37.9 and (−30.7, −54.1), respectively, and total < 0 across the virtual population. Given the initial conditions (15), the latter follows because SL tumours fully occupy the available space at the start of treatment. Tumours in the NL cohort typically respond less well to treatment, with larger median and (Q1, Q3) values of viable equal to −4.57 and (−3.52, −6.53), respectively, and total > 0, for at least 90% of tumours. When the net RT-induced cell death is minimal, NL tumours, which do not occupy all available free space at the start of treatment, can grow larger due to increases in the dead cell volume. We also note that the value of total − viable is larger for SL tumours since they accumulate more dead material.
In both regimes, we observe outliers, which undergo much larger reductions in T + T S and than the average tumour. This suggests that certain parameter values Fig. 6 The scatter plots show the values of viable and total , following a conventional fractionation schedule, for each (q 1 , q 3 ) pair used to generate the set of virtual NL tumours. viable and total increase with q 3 and decrease with q 1 (Colour figure online) within the NL and SL regimes correspond to tumours which are more sensitive to RT than the average NL and SL tumour. The influence of the oxygen consumption rates, q 1 and q 3 , on treatment outcome following a conventional fractionation schedule. We now investigate the role of q 1 and q 3 in tumour response to RT. The scatter plots in Fig. 6 show the values of viable and total across the (q 1 , q 3 ) pairs which define the NL virtual population. The response of NL tumours is most sensitive to the value of q 3 , with smaller values leading to greater reductions in viable and total cell volumes. Further, higher values of q 1 are also associated with larger reductions in viable and total cell volumes. To understand these findings, we study the response to RT of four representative tumours corresponding to (q 1 , q 3 ) sets, A 1 , B 1 , C 1 and D 1 (see Fig. 6 and Table 2).
Comparing the response of tumours A 1 and B 1 to tumours C 1 and D 1 in Fig. 7a, we see that a smaller value of q 3 implies higher average oxygen levels and slower cell proliferation (since q 2 = 0.01q 3 ). We conclude that two mechanisms could explain the increased efficacy of RT for low values of q 3 : (i) higher rates of RT cell kill due to increased oxygenation or (ii) limited regrowth between RT fractions due to decreased proliferation.
While oxygen levels are higher in tumours A 1 and B 1 than tumours C 1 and D 1 , respectively, when RT is applied (see Fig. 7a), their values of %T R are slightly smaller (16.65% (A 1 ) vs. 18.42% (C 1 ) and 13.79% (B 1 ) vs. 16.29% (D 1 ); see Fig. 7b). This suggests that the net increase in oxygen levels when values of q 3 are small does not significantly impact the rates of cell kill due to RT. By contrast, Fig. 7a shows that the relative increase in the viable cell volume of tumours A 1 and B 1 between fractions is marginal, whereas the relative increase in the viable cell volume of tumours C 1 and D 1 between fractions is sufficiently large for it to return to its initial volume over the weekend break from RT. This indicates that the value of q 3 impacts the reduction in the tumour burden by modulating tumour regrowth between fractions (rather than by increasing RT-induced cell death). Figure 7b also shows that a larger value of q 1 can slightly increase the magnitude of the reductions in viable and total . Since high values of q 1 lead to lower average oxygen levels (Fig. 7a), RT cell kill rates are smaller, while the rate of cell death due to hypoxia is larger than for low values of q 1 . The balance between these two processes  17). Tumours with small values of q 3 (A 1 , B 1 ) undergo sustained reductions in T + T S during treatment, leading to larger reductions in total volume. Tumours with large values of q 3 (C 1 , D 1 ) experience transient reductions in T + T S as they regrow to their initial volume during the weekend break from RT, and their total volume increases as they accumulate dead material. A high value of q 1 (B 1 , D 1 ) leads to a modest improvement in the tumour response, but does not yield a large reduction in tumour volume (Colour figure online) determines whether cell death increases or decreases as q 1 increases. For tumours C 1 and D 1 , Fig.7a shows that the reduction in T + T S following RT is greater and the increase in T R is smaller for larger values of q 1 . This confirms that a larger reduction in tumour burden can be achieved for large values of q 1 despite a reduction in RT-induced cell death: in such cases, increased cell death due to hypoxia drives the reduction in tumour volume. Overall, we have shown that both low values of q 3 and high values of q 1 characterise the best NL responders. Since Fig. 7a, b suggest that the value of q 1 has a less significant influence on tumour reduction than q 3 , we conclude that growth limitation between RT fractions, rather than high rates of cell death due to RT or oxygen insufficiency, has the greatest influence on the efficacy of RT for NL tumours.   Figure 8 shows the values of viable and total across the (q 1 , q 3 ) pairs which define the SL virtual population. The response of SL tumours is sensitive to the values of both q 1 and q 3 : greater reductions in viable cell volume are obtained for smaller values of q 1 and/or q 3 , while greater reductions in total cell volume are obtained for smaller values of q 3 . To understand these results, we study the response to RT of tumours corresponding to four representative (q 1 , q 3 ) sets A 2 , B 2 , C 2 and D 2 (see Fig. 8 and Table 3).
Figures 9a, b reveal that A 2 and C 2 accumulate a larger number of dead cells than tumours B 2 and D 2 . This difference in tumour composition is amplified during treatment and the parameter which influences this distinction most is q 1 . Figure 9a shows that for low values of q 1 (tumours A 2 and C 2 ), the intratumoural oxygen concentration, c, is at least 10-fold higher than for high values of q 1 (tumours B 2 and D 2 ). In particular, c c min throughout treatment when q 1 is small, which means that there is no cell death due to nutrient insufficiency and cell death is solely attributable to RT. Therefore, the decrease in viable cell volume (and corresponding increase in dead cell volume) in tumours A 2 and C 2 following each RT fraction is driven by cell kill due to RT, which is enhanced by low values of q 1 .
For tumours B 2 and D 2 , Fig. 9a also shows that, even though the oxygen concentration transiently drops below c min when each RT fraction is applied, there is a net increase in c throughout treatment and, in particular, the weekly average oxygen concentration remains above c min (result not shown). Therefore, we expect RT cell kill to increase during the fractionation schedule and cell death due to hypoxia to decrease. Since RT cell kill remains limited by low oxygen levels for both tumours, neither of the Fig. 9 (a) For a conventional fractionation schedule, we numerically solve Eqs. (10)-(13) for t ∈ (0, 8×10 4 ] subject to the initial conditions (15). In A 2 -D 2 , we fix V 0 = 0.005 and (q 1 , q 3 ) as indicated by the points A 2 , B 2 , C 2 and D 2 in Fig. 8, which correspond to SL tumours. (b) Bar graph showing the mean composition of tumours A 2 -D 2 in the last week of a conventional fractionation schedule, where %T , %T S , %T R and %V 0 are defined in Eq. (17). We observe three qualitative behaviours: (i) low q 1 (A 2 , C 2 ) is associated with high oxygen levels, high rates of RT cell kill and a large accumulation of dead cell material; (ii) high q 1 and low q 3 (B 2 ) implies a smaller reduction in viable volume due to smaller RT cell kill and a greater reduction in total volume due to limited inter-fraction tumour growth and smaller dead cell accumulation; (iii) high q 1 and q 3 (D 2 ) leads to modest reductions in viable and total volumes due to low rates of RT cell kill and high rates of cell proliferation. Overall, the values of q 1 and q 3 influence the tumour composition and the total tumour volume, respectively (Colour figure online) two proposed cell death mechanisms is responsible for the increased RT efficacy for tumour B 2 compared to tumour D 2 . However, Fig. 9a reveals that the relative increase in T + T S between fractions is smaller for tumour B 2 , which is characterised by low q 3 . We, therefore, conclude that the increased RT efficacy is driven by reduced tumour regrowth between fractions (similarly to NL tumours with low q 3 ). Figure 9b further shows how low values of q 3 enable greater reductions in total tumour volume, . Since, for tumours B 2 and D 2 , the values of %T R are comparable while the value of %T is smaller for tumour B 2 , the larger reduction in observed for tumour B 2 is due to increased net cell death (as described above). In contrast, for tumours A 2 and C 2 , the values of %T are comparable while the value of %T R is smaller for tumour A 2 . The larger reduction in observed for tumour A 2 is, therefore, due to a smaller accumulation of dead material, which occurs when lower viable cell volumes (caused by slower tumour regrowth between fractions) and/or lower oxygen levels reduce RT-induced cell death.
Overall, we have shown that two mechanisms can contribute to the increased efficacy of RT for certain tumours in a SL regime. These mechanisms are cell death due to RT and limited tumour regrowth between RT doses. Their relative contributions depend on the values of q 1 and q 3 . More specifically, when q 1 is small, RT cell kill is the dominant mechanism contributing to increased net cell death and, when q 3 is also small, limited regrowth between fractions ensures a larger reduction in total tumour volume. When q 1 is large and q 3 is small, limited regrowth between fractions determines the response to RT by ensuring larger reductions in viable and total cell volumes. The effect of the dosing schedule on typical tumour response. We now consider how, for a fixed total dose, the dose rate, R, and the number of fractions per week, N f rac , affect tumour response to RT. For the virtual cohorts of NL and SL tumours, Figs. 10 and 11, respectively, show the distributions of viable and total for fractionation schedules with R ∈ {0.1, 0.2, 0.3, 0.4, 0.5} and N f rac ∈ {1, 3, 5}. For NL tumours, the mean reduction in viable volume and the difference between the viable and total volumes increase with R and N f rac . However, the maximum reduction in viable and total volumes typically decreases with R (for fixed N f rac ), and the mean and maximum increases in total volumes also increase with R and N f rac . Therefore, a higher dosing frequency and/or dose may not lead to greater RT efficacy for tumours in the NL regime. By contrast, for SL tumours, on average, the reduction in the viable and total volumes and the difference between the viable and total volumes increase with R and N f rac . The response of SL tumours is, thus, consistent with the current, standard approach to RT protocol design, which aims to maximise RT cell kill by applying the highest tolerable total dose, in sufficiently frequent fractions, to the tumour. This result is also supported by other modelling approaches, e.g., Lewin et al. (2018) developed a spatially resolved model of avascular tumour growth and RT cell death which predicted that there is a minimum RT dose, for a fixed dosing frequency, and a minimum dosing frequency, for a fixed RT dose, below which tumours grow during treatment.

Tumours in the Bistable Regime
Typical response to a conventional fractionation schedule. Figure 12a shows the average response of a tumour in a BS regime to a conventional fractionation schedule. RT has a detrimental effect as tumour regrowth between fractions and over the weekend outweighs RT-induced cell death. The dead cell volume also increases throughout treatment, implying an increase in total volume. Figure 12b further shows that, for the  (15). We set (q 1 , q 3 , V 0 ) = (0. 787, 8.38, 0.00275). This tumour represents the typical behaviour in a BS regime. (b) Violin plots representing the distributions of viable and total . While the effect of RT is deleterious for most tumours, with several outliers experiencing larger than average increases in viable and total volumes, there are tumours that exhibit larger than average decreases in viable volume (Colour figure online) BS virtual cohort, viable > 0 for at least 80% of tumours and total > 0 for all tumours. This reveals that most tumours in the BS virtual cohort respond badly to RT.
The results in Fig. 12b also indicate that a few virtual tumours are more or less sensitive to RT than the average tumour in the BS virtual population: while their total volume increases during RT, their viable volume undergoes a 20 − 80% decrease or 40 − 80% increase, respectively, by the end of treatment. We investigate the response to RT of these outliers in more detail in the following section. The influence of q 1 , q 3 and V 0 on treatment outcome following a conventional fractionation schedule. As for tumours in monostable regimes, we study the influence of q 1 and q 3 on tumour response to RT, but we also study the role played by the vascular volume, V 0 . More specifically, we introduce a function V d , which quantifies how close where V N and V S are the threshold values of V 0 below and above which only NL and SL steady states exist, respectively. Further, In particular, V d 0 for tumours which are close to the NL regime, whereas V d 1 for tumours which are close to the SL regime. The scatter plots in Fig. 14 show the values of viable and total across the (q 1 , q 3 ) and (q 1 , V d ) pairs corresponding to the BS virtual population. We note that the values of q 3 and V d are correlated: for fixed q 1 , the lowest value of V d corresponds to the highest value of q 3 and vice versa. It is, therefore, sufficient to describe the response of tumours in a BS regime with respect to the values of q 1 and V d . The largest reductions in viable volume are obtained for lower values of q 1 and V d 1, whereas the largest increases in viable volume are obtained for higher values of q 1 and intermediate values of V d . Those tumours with the smallest and largest values of viable also undergo the largest increases in total volume: for intermediate to high values of V d , total decreases as V d and q 1 increase.
We now select four representative (q 1 , q 3 , V d ) sets A 3 , B 3 , C 3 and D 3 (see Fig. 14 and Table 4) and study the corresponding tumours' response to RT. Figures 15a, b show that tumours A 3 and B 3 decrease in viable volume, with A 3 experiencing a larger than average reduction, while C 3 and D 3 increase in viable volume, with C 3 experiencing   (19), this suggests that the behaviour of tumours in a BS regime that lie sufficiently close to the SL or NL regions will be, respectively, similar to that of SL or NL tumours with values of q 1 and q 3 of the same order of magnitude. More specifically, tumours A 3 and B 3 respond to RT similarly to SL tumours C 2 and D 2 , respectively (recall Fig. 9a, b), while tumour D 3 responds similarly to NL tumour D 1 (recall Fig. 7a, b). For tumour A 3 , this response involves an initial large increase in viable and total volume as the tumour evolves towards its SL steady state, followed by a substantial increase in RT cell kill, the average oxygen concentration and the dead cell volume. Despite the reduction in viable volume, the accumulation of dead material implies a significant increase in total volume. The same qualitative  (17). A 3 and B 3 decrease in viable volume and increase in total volume, while C 3 and D 3 increase in both viable and total volumes. Tumour A 3 , which has q 1 and V d ≈ 1, experiences larger than average decreases and increases in viable and total volumes, respectively. By contrast, tumour C 3 , which has an intermediate value of V d , experiences larger than average increases in viable and total volumes (Colour figure online) behaviour is observed for tumour B 3 , with less RT-induced cell death and dead material accumulation as the oxygen concentration remains significantly lower than for A 3 . As a result, the increase in total volume is also smaller. For tumour D 3 , cell death due to RT and hypoxia is outweighed by the tumour regrowth between fractions, leading to small increases in viable and total volumes. Further, while tumour C 3 lies closest to the SL region (V d ≈ 0.8), it does not transition from the basin of attraction of its NL steady state to its SL steady state, unlike tumours A 3 and B 3 . In particular, the increase in the oxygen concentration for C 3 is not sufficiently rapid for the tumour to enter, during treatment, a SL regime where, on average, c > c min . Therefore, the increase in viable volume is constant, but gradual, with a smaller accumulation of dead material. This explains why C 3 undergoes a larger than average increase in viable volume, with a moderate increase in total volume.
In summary, we have identified two extremal regions of parameter space in which tumours in a BS regime undergo larger decreases or increases in viable volume (and larger increases in total volume) than the typical tumour in this regime. Tumours which are sufficiently near to the boundary of the BS and SL regimes and consume little oxygen for maintenance experience larger than average decreases in viable cell volume as RT cell death is enhanced by higher oxygen levels. By contrast, tumours which are close to the boundary between the BS and SL regimes, but not sufficiently close, undergo larger than average increases in viable volume, regardless of the value The mean and maximum values of viable increase, while its minimum value decreases, as R and N f rac increase. The mean, minimum and maximum values of total increase as R and N f rac increase. There is an exception for N f rac = 5, where the maximum value of viable and minimum value of total decrease with R ≥ 3 (Colour figure online) of q 1 . This occurs as they attempt and fail to transition from their NL steady state to their SL steady state and, thus, RT cell death remains limited by low oxygen levels and outweighed by tumour regrowth between fractions. The effect of the dosing schedule on typical tumour response. For the virtual population of tumours in a BS regime, we show in Fig. 16 how the dose rate, R, and the number of fractions per week, N f rac , affect tumour response to RT when the total dose is fixed. On average, a higher number of fractions per week (for fixed R) and a higher dose rate (for fixed N f rac ) lead to greater increases in the viable and total cell populations. While these results contrast with those for tumours in SL regimes, we see that the maximum reduction in viable volume increases with R and N f rac , similarly to SL tumours. Overall, these results indicate that, in most cases, increasing the RT dose and frequency may be deleterious (similarly to NL tumours).

Post-treatment Tumour Growth Dynamics
In the previous section, we discussed the short-term response to RT of tumours in different growth regimes, distinguishing between tumours in monostable (NL and SL) and bistable regimes. We now investigate the long-term response to RT by studying post-treatment tumour growth dynamics and, in particular, the tumour steady states attained following treatment.

Steady State Analysis
We first perform a steady state analysis of the system (10)-(13) to understand the potential long-term effects of RT. Upon completion of a radiation protocol, we have R ≡ 0 thereafter. We, therefore, seek steady state solutions by setting R = 0 and d dt = 0 in Eqs. (10)-(13) and solving the following system Denoting the steady state solutions by T * , T * S , T * R and c * , respectively, Eq. (22) implies that T * S = η R ξ T * R . Therefore, we have either T * S = T * R = 0 or T * S , T * R > 0. Suppose that T * S , T * R > 0. We can show, by contradiction, that there are no physically realistic steady state solutions satisfying this condition by, first, proving that there are no SL steady states with T * S , T * R > 0 and, then, proving that there are no NL steady states with T * S , T * R > 0. If c * ≥ c min , then Eq. (20) gives since T * S > 0 and * < 1 by assumption and μ > 0 and q 2 > 0 by definition. Since T * > 0 is required for a physically realistic solution, there are no SL steady states with T * S , T * R > 0. If 0 < c * < c min , then Eq. (20) implies that Since q 2 = δ 1 , we have Then, Eq. (21) implies that Since θ 2 q 2 = δ 1,S , we have Comparing Eqs. (26) and (28), we obtain a contradiction. This implies that there are no NL steady states with T * S , T * R > 0. We, therefore, conclude that NL and SL steady state solutions of the system (10)-(13) must have T * S = T * R = 0. It is then straightforward to show that the solutions of the system (20)-(23) with T * S = T * R = 0 are equal to the steady state solutions in the absence of treatment (Colson et al. 2022) (see Appendix A).
We have shown that RT preserves the steady states and growth regimes observed in the absence of treatment. We conclude that, given T (0) = T * , tumours in monostable regimes at the start of treatment will return to their original tumour volume, 0 = Fig. 17 For a conventional fractionation schedule, we numerically solve Eqs. (10)-(13) for t ∈ (0, 2.5×10 5 ] subject to the initial conditions (15). In A 3 -D 3 , we fix V 0 = 0.00275 and (q 1 , q 3 ) as indicated by the points A 3 , B 3 , C 3 and D 3 in Fig. 14, which correspond to tumours in a BS regime. Tumours C 3 and D 3 evolve to their NL steady states following treatment, whereas tumours A 3 and B 3 switch to their SL steady states (Colour figure online) T * + V 0 , and composition (T * S = T * R = 0) after RT. In contrast, tumours in a BS regime either return to the original, NL steady state, or evolve to the SL steady state.

Elucidating Conditions for RT to Drive Steady State Switching of Tumours in Bistable Regimes
The steady state analysis showed that tumours in a BS regime may attain either a NL or a SL steady state following treatment. In particular, such tumours may undergo large increases in tumour volume in response to RT as they switch from a NL steady state to a larger SL steady state. Recall the tumours A 3 -D 3 that we defined in Sect. 4.2.2: Fig. 17 shows their response to RT both during, and following, a conventional fractionation schedule.
Tumours C 3 and D 3 underwent increases in viable volume during treatment and then returned to their NL steady states following treatment: the effect of RT was not strong enough to cause a switch in steady state. By contrast, tumours A 3 and B 3  Table 4. Comparing (a), (b) and (c), we observe that lower RT doses and less frequent dosing both prevent the tumour C 3 from evolving to the SL steady state following treatment (Colour figure online) experienced reductions in viable volume during treatment and then evolved to their SL steady states following treatment. The oxygen concentration in both of these tumours increased beyond the hypoxic threshold, c min , during (A 3 ) or following (B 3 ) RT and remained above this threshold level thereafter. This enabled the viable cell population to grow unchecked until the SL equilibrium was reached.
In contrast to tumours C 3 and D 3 , we recall that tumours A 3 and B 3 are characterised by V d ≈ 1, where V d is defined in (18). They are also, respectively, characterised by high and low values of q 1 , the oxygen consumption rate for maintenance. This suggests that tumours which are near to the boundary between the BS and SL regions in parameter space are most susceptible to undergoing a switch in steady state volume in response to RT, irrespective of the value of q 1 . This observation holds across a range of RT protocols (see Appendix C).
We now consider how the dosing schedule affects the long-term dynamics of tumours in a BS regime. In Fig. 18a, b, we show the response of tumour C 3 to two fractionation protocols comprising either 2 or 3 Gy fractions applied 5 times per week for 8 or 5.2 weeks, respectively. A switch in steady state is observed for 3 Gy fractions. This suggests that the likelihood of a tumour switching steady state increases with dose, a consistent trend in our numerical study (see Appendix C). Figure 18c additionally shows the response of tumour C 3 to 3 Gy fractions applied 3 times per week for 8.67 weeks. Comparing this figure to Fig. 18b highlights how a lower dosing frequency can prevent the transition from a NL to a SL steady state for tumours in BS regimes (see Appendix C).
These results suggest that a lower RT dose and dosing frequency may prevent uncontrolled increases in tumour volume following RT for tumours in BS regimes. As with our observations for short-term treatment responses, this challenges the assumption that a higher dose, applied with a higher frequency, will lead to a greater reduction in tumour volume.

Discussion
Cancer is a heterogeneous disease. In particular, tumours can exhibit widely varying responses to treatments. As a result, the success of existing therapies, which are typically applied following a "one-size-fits-all approach", can be highly variable. Patient-specific treatment design could aid in overcoming these barriers to treatment efficacy, but this requires increased understanding of the factors which affect tumour sensitivity to treatment. In this paper, we investigated how two distinct mechanisms of growth arrest can influence tumour responses to radiotherapy (RT).
We extended an existing model of solid tumour growth which distinguishes between nutrient limited (NL) and space limited (SL) growth control (Colson et al. 2022). In the absence of treatment, this model exhibits three growth regimes: (i) NL, where a tumour attains a NL steady state at which cell proliferation and death balance, (ii) SL, where a tumour attains a SL steady state when cell proliferation ceases due to space constraints, with no cell death, and (iii) bistable (BS), where stable NL and SL steady states coexist. In this paper, we investigated how tumours in each regime respond to RT. We found that the short-and long-term responses of tumours in monostable regimes (i.e., NL and SL) can be distinguished from those of tumours in BS regimes.
Tumours in the SL regime typically respond well to RT in the short-term, as both their viable and total volumes decrease during fractionation, while tumours in the NL regime typically respond less well, since their total volume increases despite a reduction in viable volume. However, certain NL and SL tumours respond significantly better than the average tumour in their respective regimes. By identifying parameter regions which give rise to these outliers, we determined different mechanisms that underpin successful RT. For NL tumours, RT efficacy is maximised when regrowth between fractions is minimised, while, for SL tumours, increased RT efficacy may be due to limited regrowth (as for NL tumours) and/or RT cell kill. The additional SL-specific mechanism is a consequence of low rates of RT cell kill for NL tumours due to low oxygenation. This explains how the different growth arrest mechanisms that characterise the NL and SL regimes can affect short-term tumour response to RT. In the long-term, tumours in the NL and SL regimes always return to their pretreatment steady state volume, irrespective of the effects of RT during treatment. Our model, therefore, predicts that any change in the tumour burden during radiation is transient for these tumours. This result follows from our simplifying assumption that the vascular volume remains constant during tumour growth and RT. We explain below how we could extend the model to relax this assumption.
We also found that most tumours in the BS regime respond badly to RT in the shortterm, as their viable and total cell volumes increase during RT. As for monostable regimes, outliers which lie, in parameter space, near the boundary between BS and SL regions, exhibit more extreme responses to RT. In these cases, the intratumoural oxygen concentration is close to, and smaller than, c min , the threshold concentration below which cells die due to nutrient insufficiency. If RT induces a net increase in oxygen levels such that c > c min , then cell death due to nutrient insufficiency ceases and RT drives the tumour to its SL steady state. This leads to a significant increase in RT-induced cell death and dead cell accumulation, resulting in large decreases and increases in viable and total volumes, respectively. By contrast, if RT induces a net increase in oxygen levels such that c ≤ c min , then RT causes large increases in viable and total volumes as the tumour grows towards, and fails to reach, its SL steady state. Here, RT cell kill is outweighed by tumour growth between fractions throughout treatment. Irrespective of whether these outliers experience increases or decreases in viable volume, they evolve to their larger SL steady state following RT. Therefore, the model predicts that, in the BS regime, RT usually has a detrimental effect on tumour growth.
A final key result relates to RT dosing schedules. We found that, in SL regimes, applying larger doses at higher frequency typically increases RT efficacy, whereas, in NL and BS regimes, administering lower doses at lower frequency can increase RT efficacy for outliers and lessen, or prevent, large increases in tumour burden across the virtual cohorts. The latter is a counter-intuitive result and challenges the assumption that giving the maximum tolerable dose is the best course of treatment. In practice, we are unlikely to know in which growth regime a patient's tumour lies when treatment starts. It would be interesting, in future work, to investigate whether we can determine a tumour's growth regime by monitoring its response to a given treatment protocol. If we can establish that a tumour is in a SL regime, then this would allow us to adapt the treatment protocol to maximise the reduction in tumour burden, e.g., by increasing the RT dose or dosing frequency. Alternatively, if a tumour is in a NL or BS regime, then it might be preferable to halt treatment early in order to prevent large increases in tumour burden.
While some of our findings support experimentally and clinically observed phenomena, others challenge some common assumptions. The significance of these results is limited by the fact that they remain to be tested against relevant clinical and/or experimental data. In future work, we will aim to calibrate and validate our model using a range of experimental data, with the caveat that it may be difficult to estimate all of the parameters from the data. For instance, in vitro assays involving 2D monolayers, where all tumour cells have abundant nutrient and proliferate until they reach confluence (Kapałczyńska et al. 2018), could be used to test our predictions for SL tumours. Further, in vitro assays involving 3D spheroids, which reach an equilibrium at which cell proliferation and cell death rates balance (Costa et al. 2016), could be used to test our predictions for NL tumours. Finally, if we can determine a tumour's growth regime from its initial response to RT, then, as mentioned above, we could also use in vivo tumour data from mice treated with RT to test model predictions across regimes.
In addition to testing our model against experimental data, in future work we aim also to relax some of its simplifying assumptions. First, we will view the vascular volume as a dynamic variable, which evolves in response to treatment and angiogenic signalling. In doing so, we will obtain a more realistic description, not only of the co-evolution of the tumour cells and vasculature, but also of responses to treatments which affect tumour and endothelial cells. In particular, if the vascular volume evolves over time, then treatment may alter the steady state tumour volumes. Second, since we found that the response of certain tumours (e.g., SL tumours) improves with higher doses applied at higher frequency, it would be interesting to incorporate a trade-off between anti-tumour effects and detrimental side-effects. For example, we could model the damage experienced by surrounding healthy tissue following each RT dose and impose a maximum value of cumulative damage over the RT protocol, above which the RT toxicity is considered intolerable (e.g., see Hanin and Zaider 2013;Huang et al. 2015;Kuznetsov and Kolobov 2023;Stocks et al. 2017). This would enable us to assess which dosing regimens are clinically feasible and to evaluate tumour responses to RT protocols designed using the maximum tolerated dose.
Author Contributions All authors conceived and designed the study, C.C. executed the study and all authors contributed to the writing of the manuscript.
Funding CC is supported by funding from the Engineering and Physical Sciences Research Council (EPSRC).

Conflicts of interest
We declare we have no competing interests.
Code Availability MATLAB code to generate the three virtual tumour populations and to numerically solve the ODE model is available at: https://github.com/chloeacolson/InvestigatingTumourResponsesRadiotherapy.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.  U (a, b) used to sample the pairs of q 1 and q 3 values that correspond to the NL, SL and BS cohorts, respectively.q 1 andq 3 are defined by Eqs. (B8) and (B9), respectively Cohort NL SL BS Distribution used to sample q 1 U (q 1 , 10) U (0.01,q 1 ) U (0.01,q 1 ) Distribution used to sample q 3 U (0.01, 10) U (0.01,q 3 ) U (q 3 , 10) sets which define the three cohorts considered in this paper are available at: https:// github.com/chloeacolson/InvestigatingTumourResponsesRadiotherapy.

Appendix C Numerical Results: Steady State Switching of Tumours in the Bistable Regime
For dosing regimens with R ∈ {0.1, 0.3, 0.5} and N f rac ∈ {1, 3, 5}, the scatter plots in Fig. 19 highlight (in red) the (q 1 , V d ) pairs which correspond to tumours in a bistable regime that switch steady state. We observe that tumours which switch steady state typically have larger values of V d . We note also that the number of tumours which switch steady state increases with R and N f rac . These results are consistent with those presented in Sect. 4.3.2 for tumour C 3 .

Fig. 19
For the virtual BS tumour population and fractionation schedules with R ∈ {0.1, 0.3, 0.5} and N f rac ∈ {1, 3, 5}, the scatter plots show the (q 1 , V d ) pairs that correspond to tumours that switch (red) and do not switch (black) steady state. The former are typically characterised by larger values of V d (Colour figure online)