Investigating the influence of growth arrest mechanisms on tumour responses to radiotherapy

TL;DR

This study models how different growth arrest mechanisms in tumors influence their responses to radiotherapy, revealing that tumors limited by space respond best, while bistable tumors respond worst, guiding personalized treatment strategies.

Contribution

It extends an existing tumor growth model to include two growth arrest mechanisms and analyzes their impact on radiotherapy effectiveness across different regimes.

Findings

Tumors in the space-limited regime respond best to radiotherapy.

Bistable regime tumors respond worst to radiotherapy.

Optimal dosing regimens vary depending on the growth arrest mechanism.

Abstract

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 (NL), 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.

Tables

Table 1. Table 1 : List of dimensionless parameters and their default values.

Parameters	Definition	Value(s)
$c_{min}$	Anoxic oxygen threshold	${10}^{-2}$
$g$	Rate of oxygen exchange per unit vascular volume	$5$
$k$	Conversation factor	${10}^{-2}$
$q_1$	$O_2$ consumption rate for maintenance	$[{10}^{-2},10]$
${\theta }_1$	Proportionality constant relating $q_1$ and $q_{1,S}$	$10$
$q_3$	$O_2$ consumption rate for proliferation	$[{10}^{-2},10]$
$q_2$	Proliferation rate	$q_3/k$
${\theta }_2$	Proportionality constant relating $q_2$ and $q_{2,S}$ ($q_3$ and $q_{3,S}$)	$0.1$
${\delta }_1,{\delta }_{1,S}$	Rates of death due to nutrient insufficiency	$q_2,{\theta }_2{q}_2$
$V_0$	Vascular volume	$(0,5\times {10}^{-3}]$
$\nu$	RT sub-lethal damage rate	$10$
$\lambda ,{\lambda }_S$	RT lethal damage rate	$1$
$\mu$	Repair rate constant	$5\times {10}^{-3}$
$\zeta$	Rate of death by mitotic catastrophe	$5\times {10}^{-4}$
${\eta }_R$	Clearance rate of cells killed by RT	$5\times {10}^{-5}$

Table 2. Table 2: Parameter sets A 1 subscript 𝐴 1 A_{1} , B 1 subscript 𝐵 1 B_{1} , C 1 subscript 𝐶 1 C_{1} and D 1 subscript 𝐷 1 D_{1} corresponding to the representative NL tumours.

Tumour	$q_1$	$q_3$
$A_1$	$8.91\times {10}^{-1}$	$1.14\times {10}^{-1}$
$B_1$	$7.78$	$4.01\times {10}^{-2}$
$C_1$	$8.91\times {10}^{-1}$	$9.75$
$D_1$	$7.60$	$9.94$

Table 3. Table 3: Parameter sets A 2 subscript 𝐴 2 A_{2} , B 2 subscript 𝐵 2 B_{2} , C 2 subscript 𝐶 2 C_{2} and D 2 subscript 𝐷 2 D_{2} corresponding to the representative SL tumours.

Tumour	$q_1$	$q_3$
$A_2$	$1.51\times {10}^{-1}$	$2.10\times {10}^{-1}$
$B_2$	$2.14$	$1.43\times {10}^{-1}$
$C_2$	$3.21\times {10}^{-2}$	$9.53$
$D_2$	$2.14$	$7.61$

Table 4. Table 4: Parameter sets A 3 subscript 𝐴 3 A_{3} , B 3 subscript 𝐵 3 B_{3} , C 3 subscript 𝐶 3 C_{3} and D 3 subscript 𝐷 3 D_{3} corresponding to the representative tumours in the BS cohort.

Tumour	$q_1$	$q_3$	$V_d$
$A_3$	$4.55\times {10}^{-2}$	$7.94$	$0.993$
$B_3$	$1.06$	$4.92$	$0.991$
$C_3$	$3.92\times {10}^{-2}$	$9.96$	$0.791$
$D_3$	$1.36$	$9.93$	$0.00376$

Equations

$$\text{diff}Tt=\underset{\text{shortstackrate of tumour cell proliferation}}{\underbrace{q_2^{\ast }cT(S_{max}-(T+V_0))}}-\underset{\text{shortstackrate of cell death due to nutrient starvation}}{\underbrace{\left[{\delta }_1^{\ast }(c_{min}^{\ast }-c)\right]H(c_{min}^{\ast }-c)T}},$$

$$\text{diff}ct=\underset{\text{shortstackrate of nutrient delivery}}{\underbrace{g^{\ast }(c_{max}^{\ast }-c)V_0}}-\underset{\text{shortstackbaseline rate of nutrient consumption}}{\underbrace{q_1^{\ast }cT}}-\underset{\text{shortstackadditional rate of nutrient consumption for proliferation}}{\underbrace{q_3^{\ast }cT(S_{max}-(T+V_0))}},$$

$$H(x)=\{\begin{array}{l}1,\text{if }x\ge 0,\\ 0,\text{if }x<0.\end{array}$$

$$S_{LQ}(D)=\exp \left(-\left(\alpha D+\beta D^2\right)\right),$$

$$R(t)=\{\begin{array}{l}D/{\delta }_R,\text{if}{t}_R\le t\le t_R+{\delta }_R,\\ 0,\text{otherwise.}\end{array}$$

$$\begin{gathered} \text{diff}Tt=q_2^{\ast }cT(S_{max}-\Sigma )-{\delta }_1^{\ast }(c_{min}^{\ast }-c)H(c_{min}^{\ast }-c)T \\ -\underset{\text{shortstackrate of direct lethal damage}}{\underbrace{{\lambda }^{\ast }cRT}}-\underset{\text{shortstackrate of sub-lethal damage}}{\underbrace{{\nu }^{\ast }cRT}}+\underset{\text{shortstackrate of repair}}{\underbrace{{\mu }^{\ast }{T}_S}}, \end{gathered}$$

$$\begin{gathered} \text{diff}{T}_{S}t=q_{2,S}^{\ast }c{T}_S(S_{max}-\Sigma )-{\delta }_{1,S}^{\ast }(c_{min}^{\ast }-c)H(c_{min}^{\ast }-c)T_S \\ +{\nu }^{\ast }cRT-{\mu }^{\ast }{T}_S-\underset{\text{shortstackrate of post-RT death due to MC}}{\underbrace{{\xi }^{\ast }{T}_S}}-\underset{\text{shortstackrate of indirect lethal damage}}{\underbrace{{\lambda }_S^{\ast }cR{T}_S}}, \end{gathered}$$

$$\text{diff}{T}_{R}t={\lambda }^{\ast }cRT+({\xi }^{\ast }+{\lambda }_S^{\ast }cR)T_S-\underset{\text{shortstackrate of clearance}}{\underbrace{{\eta }_R^{\ast }{T}_R}},$$

$$\text{diff}ct=g^{\ast }(c_{max}-c)V_0-q_1^{\ast }cT-q_3^{\ast }cT(S_{max}-\Sigma )-q_{1,S}^{\ast }c{T}_S-q_{3,S}^{\ast }c{T}_S(S_{max}-\Sigma ),$$

$$\hat{T}=\frac{T}{S_{max}},{\hat{T}}_S=\frac{T_S}{S_{max}},{\hat{T}}_R=\frac{T_R}{S_{max}},\widehat{V_0}=\frac{V_0}{S_{max}},$$

$$\hat{c}=\frac{c}{c_{max}},\hat{R}=\frac{R}{R_{max}},\hat{t}=\frac{t}{\tau }.$$

$$\text{diff}Tt=q_{2}cT(1-\Sigma )-({\delta }_1(c_{min}-c)H(c_{min}-c)+\lambda cR+\nu cR)T+\mu T_S,$$

$$\text{diff}{T}_{S}t={\theta }_2{q}_{2}c{T}_S(1-\Sigma )-({\delta }_{1,S}(c_{min}-c)H(c_{min}-c)+{\lambda }_{S}cR+\mu +\xi )T_S+\nu cRT,$$

$$\text{diff}{T}_{R}t=\lambda cRT+(\xi +{\lambda }_{S}cR)T_S-{\eta }_R{T}_R,$$

$$\text{diff}ct=g(1-c)V_0-q_1(T+{\theta }_1{T}_S)c-q_3\left(T+{\theta }_2{T}_S\right)c(1-\Sigma ),$$

$$q_1=q_1^{\ast }{S}_{max}\tau ,q_3=q_3^{\ast }{S}_{max}^2\tau ,q_2=q_2^{\ast }{S}_{max}^2{c}_{max}^{\ast }\tau ,k=\frac{c_{max}^{\ast }}{S_{max}}{k}^{\ast },$$

$$c_{min}=\frac{c_{min}^{\ast }}{c_{max}^{\ast }},{\delta }_1={\delta }_1^{\ast }{c}_{max}\tau ,{\delta }_{1,S}={\delta }_{1,S}^{\ast }{c}_{max}\tau ,g=g^{\ast }{S}_{max}\tau ,$$

$$\lambda ={\lambda }^{\ast }{c}_{max}{R}_{max}\tau ,{\lambda }_S={\lambda }_S^{\ast }{c}_{max}{R}_{max}\tau ,\nu ={\nu }^{\ast }{c}_{max}{R}_{max}\tau ,$$

$$\mu ={\mu }^{\ast }\tau ,\xi ={\xi }^{\ast }\tau ,{\eta }_R={\eta }_R^{\ast }\tau .$$

$$(T,T_S,T_R,c)=(T^{\ast },0,0,c^{\ast }),$$

$${\Delta }_{\text{viable}}:=100\times \frac{(\bar{T}+\overline{T_S})-T_0}{T_0}\text{and}{\Delta }_{\text{total}}:=100\times \frac{\bar{\Sigma }-{\Sigma }_0}{{\Sigma }_0},$$

$${\Sigma }_0=T_0+V_0,\Sigma =T+T_S+T_R+V_0.$$

$$\%T:=100\times \frac{\bar{T}}{{\Sigma }_0},\%T_S:=100\times \frac{\overline{T_S}}{{\Sigma }_0},\%T_R:=100\times \frac{\overline{T_R}}{{\Sigma }_0}\text{and}\%V_0:=100\times \frac{V_0}{{\Sigma }_0}.$$

$$V_d(V_0)=\frac{V_0-V_N}{V_S-V_N},$$

$$\{\begin{array}{l}{V}_d\to 0\text{ as }{V}_0\to V_N,\\ V_d\to 1\text{ as }{V}_0\to V_S.\end{array}$$

$$q_{2}cT(1-\Sigma )-{\delta }_1(c_{min}-c)H(c_{min}-c)T+\mu T_S=0,$$

$${\theta }_2{q}_{2}c{T}_S(1-\Sigma )-({\delta }_{1,S}(c_{min}-c)H(c_{min}-c)+\mu +\xi ))T_S=0,$$

$$\xi T_S-{\eta }_R{T}_R=0,$$

$$g(1-c)V_0-q_1(T+{\theta }_1{T}_S)c-q_3\left(T+{\theta }_2{T}_S\right)c(1-\Sigma )=0.$$

$$T^{\ast }=\frac{-\mu T_S^{\ast }}{q_2{c}^{\ast }(1-{\Sigma }_R^{\ast })}<0,$$

$$\frac{q_2}{{\delta }_1}(1-{\Sigma }^{\ast })-\frac{(c_{min}-c^{\ast })}{c^{\ast }}=-\frac{\mu T_S^{\ast }}{{\delta }_1{c}^{\ast }{T}^{\ast }}<0.$$

$$(1-{\Sigma }^{\ast })-\frac{(c_{min}-c^{\ast })}{c^{\ast }}=-\frac{\mu T_S^{\ast }}{{\delta }_1{c}^{\ast }{T}^{\ast }}<0.$$

Full text

1 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 [28, 29] 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 [11, 23, 27, 39] predict growth arrest due to a cessation of proliferation (with no explicit cell death), while those proposed by [10, 17, 18, 22] predict growth arrest due to the balance of cell proliferation and cell death.

In previous work [7], we developed a model of 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 more detail, we viewed the vascular volume as a parameter which influences nutrient and space availability and, therefore, the tumour’s carrying capacity. 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:

[TABLE]

[TABLE]

where

[TABLE]

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}^{*}\,($\mathrm{kg}$^{-1}$\mathrm{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 $\delta_{1}^{*}\,(${\mathrm{m}}^{3}$$\mathrm{kg}$^{-1}$\mathrm{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^{*}$ 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}^{*}\,($\mathrm{m}$^{-3}$\mathrm{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\,($\mathrm{m}$^{-6}$\mathrm{min}$^{-1})$.

Since the model (1.1)-(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 Equations (1.1)-(1.2) to account for the biological effects of RT.

Radiobiology.

RT is used to treat more than $50\%$ of cancer patients [24]. 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\,\mathrm{Gy}$ doses delivered Monday to Friday for up to $7$ weeks [1], the dose, dosing frequency and treatment duration can be varied to deliver a fixed total dose. The latter is often termed as the Maximum Tolerated Dose (MTD), i.e., the highest dose which does not cause adverse side-effects [14].

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 [16]. Intratumoural oxygen levels are a key factor influencing tumour radio-sensitivity, another being 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 [25] and vascular damage and necrosis observed at high doses [35, 37]. 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 [13, 19, 21, 32, 33], designing patient-specific RT dosing schedules [3, 4] and investigating the influence of inter- and/or intra- tumour heterogeneity on tumour sensitivity to RT [2, 6, 12, 31, 38].

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 [26]. The LQ model states that the fraction, $S_{LQ}$, of (tumour) cells that survive exposure to a dose $D$ ($\mathrm{Gy}$) of radiation is given by

[TABLE]

where $\alpha\geq 0$ and $\beta\geq 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 2D monolayer or 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 [8, 15, 30, 36]. 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. [30] and adopt a time-dependent description of RT.

Paper structure.

This paper is structured as follows. In Section 2, we extend the tumour growth model defined by Equations (1.1)-(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 Section 3. Then, we investigate the response of tumours characterised by different growth-limiting mechanisms in Section 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 Section 5, where we discuss our findings and outline possible avenues for future work.

2 Model development

2.1 The mathematical model

In this section, we incorporate tumour response to RT into the growth model (1.1)-(1.2). We follow the approach outlined in Neira et al. [30] 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 sub-lethally and lethally damaged by RT. We suppose that the tumour is exposed to a total dose $D\,(\mathrm{Gy})$ of RT at a constant rate $R$ over the time period $t_{R}\leq t\leq t_{R}+\delta_{R}\,($\min$)$ so that

[TABLE]

Let $\Sigma=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 Figure 1):

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $H$ is the Heaviside function defined in (1.3).

We assume that undamaged tumour cells, $T$, proliferate, die due to nutrient insufficiency and consume oxygen for proliferation and maintenance as in Equations (1.1)-(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 $\nu^{*}>0$ and $\lambda^{*}>0$, respectively. We further suppose that sub-lethal damage is either repaired at a constant rate $\mu^{*}>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 $\lambda^{*}_{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 $\xi^{*}>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}=\theta_{2}q^{*}_{2}>0$, with $\theta_{2}\in(0,1)$. The latter ensures that damaged cells proliferate more slowly than undamaged cells as they expend more energy repairing RT damage than proliferating. Accordingly, sub-lethally damaged cells consume oxygen for maintenance at a rate proportional to $c$, with rate constant $q^{*}_{1,S}=\theta_{1}q^{*}_{1}>0$, where $\theta_{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 $k>0$ such that $q^{*}_{3,S}=\frac{q^{*}_{2,S}}{k}$. Here, we assume the same conversion factor for $T$ and $T_{S}$ cells, for simplicity. Since $q^{*}_{2,S}=\theta_{2}q^{*}_{2}$ and $q^{*}_{3}=\frac{q^{*}_{2}}{k}$, we also have $q^{*}_{3,S}=\theta_{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 $\delta^{*}_{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 $\eta^{*}_{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.

2.2 Non-dimensionalisation

We non-dimensionalise Equations (2.1)-(2.5) by introducing the following scalings:

[TABLE]

The timescale of interest is fixed to be $\tau=1\,$\min$$, which is the timescale for the duration of RT, and the maximum dose rate $R_{max}=1\,\mathrm{Gy}/\mathrm{min}$ [20]. Then, given that $q^{*}_{1,S}=\theta_{1}q^{*}_{1}$, $q^{*}_{2,S}=\theta_{2}q^{*}_{2}$ and $q^{*}_{3,S}=\theta_{2}q^{*}_{3}$ and dropping hats for notational convenience, we obtain the following dimensionless system:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where we have introduced the following dimensionless parameter groupings:

[TABLE]

2.3 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 $\nu$, $\lambda$, $\lambda_{S}$, $\mu$ and $\eta_{R}$ are based on values found in the literature [30, 34] and we assume $\lambda=\lambda_{S}$, for simplicity. We also set $\xi=5\times 10^{-4}$ so that cells that undergo mitotic catastrophe have a half-life of approximately $24\mathrm{h}$. Here, we implicitly assume that the average duration of the cell cycle in healthy cells [5] and cancer cells are approximately the same. The parameters that define the RT dosing schedules (e.g. the dose rate, $R$) are summarised in the Methods section.

Further, we define tumour growth parameters as in our previous work [7], with the additional simplifying assumption that, as for the undamaged cells, $q_{1,S}=q_{2,S}$. We also set $\theta_{1}=10$ and $\theta_{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.

3 Review of the key model dynamics in the absence of treatment

In this section, we summarise the model behaviour in the absence of treatment. Setting $R\equiv 0$ in Equations (2.6)-(2.9), we recover the dimensionless form of Equations (1.1)-(1.2). In [7], we showed that this model admits two non-trivial, stable steady states (see Appendix A):

a nutrient limited (NL) steady state, attained when cell proliferation balances cell death due to nutrient starvation; 2. 2.

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\leq T<1-V_{0}$ and either $0\leq c<c_{min}$, for the NL steady state, or $c\geq 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. For fixed values of $c_{min}=0.01$, $g=5$ and $k=0.01$, Figure 2 shows these regions in $(q_{3},V_{0})$-space for three values of $q_{1}$.

Given $q_{1}$, there exists a threshold value of $V_{0}$, say $V_{N}$, which is independent of $q_{3}$, such that only the NL steady state exists for $0<V_{0}\leq V_{N}$. Tumours in this region of parameter space, e.g. $A_{1}$ and $A_{2}$, are said to be in an NL growth regime. Further, given $q_{1}$ and $q_{3}$ sufficiently large relative to $q_{1}$, there exists another threshold value of $V_{0}$, say $V_{S}$, such that only the SL steady state exists for $V_{0}\geq V_{S}>V_{N}$.Tumours in this region of parameter space, e.g. $C_{1}$, $B_{2}$ and $C_{2}$ are said to be in an SL growth regime. In this case, for $V_{N}<V_{0}<V_{S}$, the NL and SL steady states co-exist. Thus, in this region of parameter space, a tumour, 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}\lesssim 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^{*}$) 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)\ll T^{*}$, which we use to simulate tumour growth. As the values of $T^{*}$ and $c^{*}$ for NL tumours increase with $V_{0}$, tumours in a BS regime will grow to larger volumes than tumours in an 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\equiv 0$), Equations (2.6)-(2.9) describe three possible scenarios for tumour growth: (i) nutrient limited growth, where the tumour grows to an 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)\ll T^{*}$). In Section 4, we investigate how tumours in these growth regimes respond to RT in order to understand how different growth arrest mechanisms may influence tumour response.

4 Investigating tumour response to radiotherapy

4.1 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, which correspond to the NL, SL and BS regimes, respectively.

We then define the RT protocols of interest. We vary the dose amount $D\in\llbracket 0,5\rrbracket\,($\mathrm{G}\mathrm{y}$)$ and the number of doses per week $N_{frac}=\{1,3,5\}$. We assume that each dose is administered in $\delta_{R}=10\,$\mathrm{m}\mathrm{i}\mathrm{n} and, therefore, we vary the dimensionless dose rate, $R:=\frac{D}{\delta_{R}R_{max}}\in\llbracket 0.1,0.5\rrbracket$. 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\,$\mathrm{G}\mathrm{y} (or the closest multiple of $D$ to $80\,$\mathrm{G}\mathrm{y}$$).

For each set of tumours and each RT protocol, we solve Equations (2.6)-(2.9) numerically for $t\in(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 [9]. For simplicity, we impose the initial conditions

[TABLE]

where $T^{*}$ and $c^{*}$ are 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 record $\bar{T},\,\bar{T_{S}}$ and $\bar{T_{R}}$, the mean undamaged, damaged and dead cell volumes in the last week of treatment. We also define the percent change in (mean) viable and total cell volumes between the start and the end of treatment as follows

[TABLE]

where

[TABLE]

We also quantify the end-of-treatment tumour composition (relative to the total tumour volume at the start of treatment) as follows

[TABLE]

We note that the variables defined in (4.3) can be used to describe $\Delta_{total}=(\%T+\%T_{S}+\%T_{R}+\%V_{0})-100$. Finally, we record $\bar{c}$, the mean oxygen concentration in the last week of treatment, and the post-treatment steady state values of all the dependent variables.

4.2 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\times 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.

4.2.1 Tumours in monostable regimes: the NL and SL virtual tumour populations

Typical responses to a conventional fractionation schedule.

Figure 4 shows the response of two NL and SL tumours to RT, the viable tumour cell volume, $T+T_{S}$, of both decreasing 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 $\Delta_{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, and hence there is a higher rate of RT cell kill and greater accumulation of dead material, $T_{R}$.

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 $\Delta_{viable}$ and $\Delta_{total}$, following a conventional fractionation schedule, across the NL and SL virtual populations. We note that the behaviour shown in Figure 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 $\Delta_{viable}$ equal to $-37.9$ and $(-30.7,-54.1)$, respectively, and $\Delta_{total}<0$ across the virtual population. Given the initial conditions (4.1), 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 $\Delta_{viable}$ equal to $-4.57$ and $(-3.52,-6.53)$, respectively, and $\Delta_{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 increases in the dead cell volume. We also note that the value of $\Delta_{total}-\Delta_{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 $\Sigma$ than the average tumour. This suggests that certain parameter values 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 Figure 6 show the values of $\Delta_{viable}$ and $\Delta_{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 Figure 6 and Table 2).

Comparing the response of tumours $A_{1}$ and $B_{1}$ to tumours $C_{1}$ and $D_{1}$ in Figure 7, 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 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}$ at the time of each dose of RT (see Figure 7), their values of $\%T_{R}$ are slightly smaller (see Figure 8). This suggests that the net increase in oxygen levels when values of $q_{3}$ are small does not significantly impact the proportion of cell kill due to RT. By contrast, Figure 7 shows that the viable cell volume of tumours $A_{1}$ and $B_{1}$ increases marginally ($A_{1}$) or remains approximately constant ($B_{1}$) between fractions, whereas the viable cell volume of tumours $C_{1}$ and $D_{1}$ increases significantly between fractions, returning to its initial volume over the week-end 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 8 also shows that a larger value of $q_{1}$ can slightly increase the magnitude of the reductions in $\Delta_{viable}$ and $\Delta_{total}$. Since high values of $q_{1}$ lead to lower average oxygen levels (Figure 7), 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 determines whether cell death increases or decreases as $q_{1}$ increases. For tumours $C_{1}$ and $D_{1}$, Figure 7 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 Figures 7 and 8 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 9 shows the values of $\Delta_{viable}$ and $\Delta_{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 Figure 9 and Table 3).

Figures 10 and 11 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 most this distinction is $q_{1}$. Figure 10 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\gg 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}$, Figure 10 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 two proposed cell death mechanisms is responsible for the increased RT efficacy for tumour $B_{2}$ compared to tumour $D_{2}$. However, Figure 10 reveals that $T+T_{S}$ increases at a slower rate between fractions 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 11 further shows how low values of $q_{3}$ enable greater reductions in total tumour volume, $\Sigma$. 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 $\Sigma$ 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 $\Sigma$ 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_{frac}$, affect tumour responses to RT. For the virtual cohorts of NL and SL tumours, Figures 12 and 13, respectively, show the distributions of $\Delta_{viable}$ and $\Delta_{total}$ for fractionation schedules with $R\in\llbracket 0.1,0.5\rrbracket$ and $N_{frac}\in\{1,3,5\}$. For SL tumours, we see that, on average, the reductions in the viable and total volumes and the difference between the viable and total volumes increase with $R$ and $N_{frac}$. 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 a 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. [21] 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. For NL tumours, the mean reduction in viable volume and the difference between the viable and total volumes also increase with $R$ and $N_{frac}$. However, the maximum reductions in viable and total volumes typically decrease with $R$ (for fixed $N_{frac}$), and the mean and maximum total volumes also increase with $R$ and $N_{frac}$. Therefore, a higher dosing frequency and/or dose may not lead to greater RT efficacy.

4.2.2 Tumours in the bistable regime

Typical response to a conventional fractionation schedule.

Figure 14(a) 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 week-end outweighs RT-induced cell death. The dead cell volume also increases throughout treatment, implying an increase in total volume. Figure 14(b) further shows that, for the BS virtual cohort, $\Delta_{viable}>0$ for at least $80\%$ of tumours and $\Delta_{total}>0$ for all tumours. This reveals that most tumours in the BS virtual cohort respond badly to RT.

The results in Figure 14(b) 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 a tumour in the BS regime lies to the NL and SL regimes (see the schematic in Figure 15):

[TABLE]

where $V_{N}$ and $V_{S}$ are the threshold values of $V_{0}$ below and above which only NL and SL steady states exist. Further,

[TABLE]

In particular, $V_{d}\gtrapprox 0$ for tumours which are close to the NL regime, whereas $V_{d}\lessapprox 1$ for tumours which are close to the SL regime.

The scatter plots in Figure 16 show the values of $\Delta_{viable}$ and $\Delta_{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 values 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 $q_{1}$ and $V_{d}$. The largest reductions in viable volume are obtained for lower values of $q_{1}$ and $V_{d}\lesssim 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 $\Delta_{viable}$ also undergo the largest increases in total volume: for intermediate to high values of $V_{d}$, $\Delta_{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 Figure 16 and Table 4) and study the corresponding tumours’ responses to RT. Figures 17 and 18 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 a larger than average increase. While tumours $A_{3}$ and $C_{3}$ have low $q_{1}$, $V_{d}\approx 1$ for $A_{3}$ and $V_{d}\approx 0.8$ for $C_{3}$. Similarly, while tumours $B_{3}$ and $D_{3}$ have high $q_{1}$, $V_{d}\approx 1$ for $B_{3}$ and $V_{d}\approx 0$ for $D_{3}$. Given (4.5), 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 Figures 10 and 11), while tumour $D_{3}$ responds similarly to NL tumour $D_{1}$ (recall Figures 7 and 8). For tumour $A_{3}$, this 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 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 rapid 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}\approx 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 rapid enough 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 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 Figure 19 how the dose rate, $R$, and the number of fractions per week, $N_{frac}$, 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_{frac}$) 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_{frac}$, 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).

4.3 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.

4.3.1 Steady state analysis

We first perform a steady state analysis of the system (2.6)-(2.9) to understand the potential long-term effects of RT. Upon completion of a radiation protocol, we have $R\equiv 0$ thereafter. We, therefore, seek steady state solutions by setting $R=0$ and $\frac{\mathrm{d}}{\mathrm{d}t}=0$ in Equations (2.6)-(2.9) and solving the following system

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Denoting the steady state solutions by $T^{*}$, $T_{S}^{*}$, $T_{R}^{*}$ and $c^{*}$, respectively, Equation (4.8) implies that $T_{S}^{*}=\frac{\eta_{R}}{\xi}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^{*}\geq c_{min}$, then Equation (4.6) gives

[TABLE]

since $T_{S}^{*}>0$ and $\Sigma^{*}<1$ by assumption and $\mu>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 Equation (4.6) supplies

[TABLE]

Since $q_{2}=\delta_{1}$, we have

[TABLE]

Then, Equation (4.7) implies that

[TABLE]

Since $\theta_{2}q_{2}=\delta_{1,S}$, we have

[TABLE]

Comparing Equations (4.12) and (4.14), 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 (2.6)-(2.9) must have $T_{S}^{*}=T_{R}^{*}=0$. It is then straightforward to show that the solutions of the system (4.6)-(4.9) with $T_{S}^{*}=T_{R}^{*}=0$ are equal to the steady state solutions in the absence of treatment [7] (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, $\Sigma_{0}=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.

4.3.2 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 Section 4.2.2: Figure 20 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 state 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}$ experienced reductions in viable volume during treatment and then evolved to their SL steady state 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}\approx 1$, where $V_{d}$ is defined in (4.4). 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 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 B).

We now consider how the dosing schedule affects the long-term dynamics of tumours in a BS regime. In Figures 21(a) and 21(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 B). Figure 21(c) 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 Figure 21(b) highlights how a lower dosing frequency can prevent the transition from NL to SL steady state for tumours in BS regimes (see Appendix B).

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.

5 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 tumour growth which distinguishes between nutrient limited (NL) and space limited (SL) growth control [7]. 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 that 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 NL and SL regimes always return to their pre-treatment 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.

We also found that most tumours in a BS regime respond badly to RT in the short-term, 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}$, 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\leq c_{min}$, 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 a 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 which growth regime a patient’s tumour lies in 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, 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, it might be preferable to halt treatment early in order to prevent large increases in tumour burden.

In this paper, we studied the effects of RT on tumour cells and neglected its effects on the tumour vasculature. In particular, we viewed the vascular volume as a parameter which influences a tumour’s carrying capacity, rather than a dynamic variable. This simplifying assumption will cease to be valid at long times when effects such as angiogenesis and vascular remodelling become important. In future work, we will extend our model to relax this assumption, and obtain a more realistic description of tumour response to treatments which affect both tumour and endothelial cells.

Appendix A Steady state solutions in the absence of treatment

A steady state analysis of the system (2.6)-(2.9) with $R\equiv 0$ was performed in [7]. There exist two SL steady states given by

[TABLE]

$\mathrm{SS}_{1}$ is unstable for all combinations of parameters, while $\mathrm{SS}_{2}$ is stable in the parameter regions in which it is an admissible solution.

There are also two NL steady states given by

[TABLE]

where

[TABLE]

and

[TABLE]

with

[TABLE]

In the regions in which $\mathrm{SS}_{3}$ and $\mathrm{SS}_{4}$ exist as admissible steady state solutions, $\mathrm{SS}_{3}$ is stable, while $\mathrm{SS}_{4}$ is unstable.

Appendix B Numerical results: steady state switching of tumours in the bistable regime

For dosing regimens with $R\in\{0.1,0.3,0.5\}$ and $N_{frac}\in\{1,3,5\}$, the scatter plots in Figure 22 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 the $R$ and $N_{frac}$. These results are consistent with those presented in Section 4.3.2 for tumour $C_{3}$.