We first investigate the behaviour of a spatially homogeneous submodel. We then reintroduce tumour cell and EC movement and use λ - ω analysis to obtain insight into the spatio-temporal dynamics. We finish by considering the response of the full system to chemotherapy.
Spatially homogeneous submodel - no therapy
We obtain our spatially homogeneous submodel by assuming that no chemotherapy is applied and by neglecting spatial effects, i.e. assuming that c = c
b
= 0 and that the variables in (12)-(14) are functions only of time. We also assume that the timescales inherent in the oxygen equation are short in comparison to the timescale for tumour growth (as for instance assumed in [15]), i.e.
Under this assumption (14) is in a quasi-steady state so
where
are dimensionless parameter groupings. In what follows we omit the asterisks and tildes in (12)-(14) and (17) for notational convenience. After substituting (17) into (12), our model reduces to the following pair of ODEs for p and v:
which can be analysed using phase plane techniques. By setting the derivative in (19) equal to zero we deduce that the p-nullclines are:
and
We note that f (v) is a straight line. From (19) it follows that ∀p > 0,
< 0 if p >f (v) and
> 0 if p <f (v). Therefore, if 1 - d
p
s
β
- d
p
≤ 0, then ∀v > 0, f (v) < 0, and no physically realistic co-existence steady state exists. Instead the tumour decreases in size ∀v > 0 (since
< 0), and the system evolves to the linearly stable tumour-free steady-state solution, (p, v) = (0, 1). So that a physically realistic co-existence steady state may exist, we require d
p
< 1/(1 + s
β
), meaning that the apoptotic rate of the tumour cells is less than the maximum rate of proliferation.
Similarly, setting (20) equal to zero we deduce that the nontrivial v-nullcline is:
It is straightforward to show that g(v) is positive only if 0 ≤ v ≤ 1 (see Figure 1 where the v-nullclines are indicated by solid red lines). Therefore, for the p- and v-nullclines to intersect and a co-existence steady state to exist, the zero of f (v), f
zero
, must fulfill
as follows from (22) (if f
zero
≥ 1, the system evolves to the linearly stable tumour-free steady state (p, v) = (0, 1)). For smaller values of f
zero
(0 <f
zero
< 1), the stability of the co-existence steady state depends on the ratio δ/η0 in (23).
Case (a): δ/η0 ≤ 1
If δ/η0 = 1, then g(v) is a straight line, g(v) = 1 - v. As shown in Figure 1 (panel A) if δ/η0 < 1, then g(v) is a monotonically decreasing function for v > 1 - δ/η0, reaching zero at v = 1 (with g(v) < 0 for v < 1 - δ/η0). Since
> 0 for p <g(v) and
< 0 for p <g(v) (as follows from (20)), we have that at the steady state V
p
< 0 and V
v
< 0, the subscripts denoting partial derivatives and the function V given by (20). Because at the steady state the p-nullcline is such that P
p
< 0 <P
v
, it follows that tr(A) = P
p
+ V
v
< 0 and det(A) = P
p
V
v
- P
v
V
p
> 0, where A = (a
ij
) is the Jacobian matrix, implying that the co-existence steady state is always linearly stable. In Figure 1 (panel A) for a case where the steady state is a stable spiral, we show the phase plane, the nullclines and a trajectory with initial conditions corresponding to the introduction of a small number of tumour cells into a tissue with vasculature at carrying capacity ((p(0), v(0)) = (0.01, 1)). The phase planes were calculated numerically using MATLAB, and equations (19)-(20) were integrated using a Runge-Kutta method.
Case (b): δ/η0 > 1
The function g(v) is concave for v ≥ 0 with g(0) = 0 and g(1) = 0 (see Figure 1, panel B). In this case oscillatory solutions are possible. In Figure 1 we present a typical phase plane for a case where the steady state is an unstable spiral and a stable limit cycle exists. The mechanism responsible for the oscillations is vessel occlusion by the tumour cells. As the tumour cells proliferate and increase in number, more vessels become occluded. Occlusion reduces the level of oxygen and hence the rate of tumour cell proliferation which results in fewer tumour cells. The reduction in tumour burden allows the vessels to recover. Vessel recovery leads to increased oxygenation and, therefore, an increase in the rate of tumour cell proliferation and the size of the tumour cell population. With the expanding tumour mass, occlusion once again becomes important, causing the cycle to repeat (see Figure 2 for a schematic explanation). In fact, the oscillations in this submodel are similar to those observed in "predator-prey" systems (see e.g. [20]); here the tumour cells act as "predators" by occluding vessels, "the prey".
We remark that it is the particular functional form of
≡ V, which appears in (20), that gives rise to oscillatory solutions; other choices of
≡ V may not admit oscillations. In particular, we require V
v
> 0 so that an unstable steady state with a limit cycle may exist. According to the Hopf bifurcation theorem, an unstable steady state with a stable limit cycle, and thus oscillations, may appear as the bifurcation parameter is varied so that the real part of a pair of complex conjugate eigenvalues of the Jacobian matrix changes sign [20]. This corresponds to a co-existence steady state changing from a stable spiral to an unstable spiral (or vice versa). Since trA = P
p
+ V
v
> 0 is a necessary condition for an unstable spiral and P
p
< 0, we need V
v
> 0 for the existence of an unstable spiral, implying that it can only exist where g(v) has positive slope (for the problem at hand, where the slope of g(v) is positive,
< 0 to the left of the nullcline and
> 0 to the right; thus V
v
> 0). In Figure 1, if the p-nullcline is moved to the right, causing the zero of f (v), f
zero
, to increase, the co-existence steady state changes stability, and the stable limit cycle disappears via a Hopf bifurcation. It is also possible to obtain a Hopf bifurcation by increasing the rate of oxygen consumption by tumour cells, σ
p
, thus decreasing the slope of the p-nullcline. We remark that the oscillatory dynamics close to the Hopf bifurcation can be approximated by the corresponding normal form of λ-ω type (see below, and section B in the Appendix (additional file 1)).
To summarise, in this section we have shown how the existence and stability of the co-existence steady state depend on the model parameters. In particular, oscillatory solutions may occur when δ, the maximum rate of occlusion, is sufficiently high; if δ exceeds the maximum rate of endothelial cell proliferation, η0, then the co-existence steady state may be unstable and periodic oscillations are possible. We remark that a high value of δ corresponds to immature vessels that are more prone to collapse.
Spatial model with random motion of tumour cells and endothelial cells - no therapy
In this section we re-introduce spatial effects (deferring the inclusion of oxygen diffusion until we investigate chemotherapy). We continue to assume that oxygen is at quasi-steady state and does not diffuse, so that s is defined in terms of v and p by (17) and equations (12)-(13) become:
In what follows we suppose that a small mass of tumour cells is introduced at one side (x = 0) of a homogeneous vascularised tissue which is otherwise devoid of tumour cells ((p(x, 0), v(x, 0)) = (0, 1) ∀x ≠ 0). This situation could describe (i) tumour cells implanted in an experiment, (ii) metastatic cells that have entered the tissue from the supporting vasculature, or (iii) cells in an avascular tumour that arose spontaneously within a tissue. Due to random movement, the tumour cells may spread to and colonise parts of the tissue that were initially unaffected. We analyse cases for which the co-existence steady state must be linearly stable (case (a) above, corresponding to δ/η0 ≤ 1), before studying cases for which it may be unstable (case (b) above, corresponding to δ/η0 > 1). In the case of an unstable steady state we note that when the random motion of endothelial cells is equal to that of tumour cells (D
v
= 1), the system lends itself to λ-ω analysis. For mathematical tractability, we carry out the λ-ω approximation for the simplified case in which the rate at which the tumour cells consume oxygen is negligible (σ
p
= 0 in (25)). Numerical simulations with D
v
≠ 1 and σ
p
≠ 0 reveal qualitatively similar dynamics.
To obtain numerical solutions to the system (25)-(26) we use the NAG routine D03PCF in Fortran [25]. This routine utilises finite differences for the spatial discretisation. Specifically, the method of lines is used to reduce the governing PDEs to a system of ODEs, which is then solved using a backward differentiation method. In each figure legend we state the mesh size Δx.
Stable co-existence steady state
When the co-existence steady state is stable, diffusion driven instabilities (DDI) may occur for certain parameter values provided that the diffusion co-efficients differ (D
v
≠ 1) [26]. Analysis of our system shows that DDI can only occur if D
v
< 1 and we have used parameter values such that the system is driven unstable. However, for brevity, we do not present those simulations or analyse this case further.
We have carried out simulations of equations (25)-(26) for several cases in which the underlying co-existence steady state is linearly stable and the random motility co-efficients for tumour and endothelial cells are equal (D
v
= 1 and DDI is not possible). In general these simulations show that when a population of tumour cells is introduced at the left-hand boundary, it invades the domain and the system evolves to its stable steady state behind the invading front. At large times the domain becomes spatially-uniform and the system reaches a stable steady state.
Unstable co-existence steady state
We now consider the case for which the spatially-uniform (no diffusion) model possesses an unstable co-existence steady state, with predator-prey-like oscillatory solutions (and with the tumour cells acting as "predators"). In [17] Sherratt studied a predator-prey reaction-diffusion model which, like ours, admits an unstable co-existence steady state, with a stable limit cycle. His analysis focussed on what happens when a small number of predators is introduced into a spatially-homogeneous population of prey. Behind the invading front, which connects the stable predator-free equilibrium ahead of the wave with the unstable co-existence equilibrium, two types of oscillatory wakes were seen: regular spatio-temporal waves and highly irregular oscillations [17]. In general, the invading front decays exponentially, with respect to the spatial variable, to the co-existence equilibrium and, if the system is close to the Hopf bifurcation, it stays there for a significant time (even though the equilibrium is unstable) before developing into a periodic wave train [27]. The stability of the resulting wave train determines whether the oscillations remain periodic (in the case of a stable wave train) or become irregular (in the case of an unstable wave train) [17, 27].
Following [17, 27] (see also [28]) we analyse equations (25)-(26) (here with σ
p
= 0 and D
v
= 1) by using normal form analysis to approximate it by one of λ-ω type near the Hopf bifurcation. The result of our analysis (for details of the analysis see section B in the Appendix (additional file 1)) is a prediction about the stability of the periodic waves behind the invading front. By keeping all other parameters fixed, and also ensuring that the system is close to Hopf bifurcation, it is possible to show that the stability of a wave train depends solely on the maximum rate of vessel occlusion, δ, and the tumour apoptosis rate, d
p
. In Figure 3, we have partitioned (δ, d
p
)-space into different regions according to whether the wave behind the invading front is stable (in the shaded region) or unstable (in the white region). We remark that once the other parameters have been fixed, the value of d
s
, the oxygen consumption rate of other cellular components, determines whether the system is near Hopf bifurcation, and thus whether the stability prediction is reliable.
In Figure 4, we present a simulation for which the λ-ω analysis, as presented in Figure 3, predicts stable periodic travelling waves behind the invading front. We note that damped oscillations connect the tumour-free steady state ahead of the invading front with the co-existence steady state, a solution which is qualitatively similar to those presented in [17]. Since the co-existence steady state is unstable, oscillations develop; in this case these are regular spatio-temporal waves because the periodic wave is stable. By comparing panels A and B in Figure 4, we see that the portion of the spatial domain in which the system is at the co-existence steady state increases over time; thus the spatial domain appears to become transiently spatially uniform. As established for a predator-prey system in [29], this occurs because the speed of the interface between the "plateau" (where the system is at steady state) and the region of regular oscillations is less than the speed of the invading front. We remark that for large times, after the invading front has reached the right-hand boundary, regular temporal oscillations persist at fixed spatial positions. After the invading front has reached the right-hand boundary it is reflected (due to the no-flux boundary conditions), and regular waves, travelling in the opposite direction to that of the original front, develop. The no-flux conditions imply that the boundaries are impermeable to the tumour cells, an assumption whose validity depends on the tissue of interest. For example, while cancer cells may be able to colonise a softer tissue, a more rigid material, such as bone, may halt cancer invasion.
In Figure 5 the waves that develop in the wake of the invading front are unstable, and irregular dynamics are observed. Since the system is close to the Hopf bifurcation, i.e. it is weakly unstable, the oscillations that develop behind the invading front are initially regular and only become irregular at later times (similar behaviour was observed in [17]). Once the cancer cells have colonised the entire tissue, irregular spatio-temporal oscillations persist.
In the above simulations oxygen diffusion was neglected and oxygen was assumed to be at quasi-steady state. The simulations were repeated for the full system, as defined by equations (12)-(14), with c = 0 in (12) and the following choice of parameter values: D
v
= 1,
and
, wherein
and
correspond to the values of d
s
and σ
p
used in (25), as follows from (18). Guided by the discussion in section A of the Appendix (additional file 1), in our simulations of the full system we fixed D
s
= 2 × 105 and h
s
= 1.6 × 106. The solutions of the full system for both stable and unstable co-existence steady states were qualitatively similar to those obtained for the simpler system (25)-(26) (and are not presented). We therefore conclude that the quasi-steady state assumption for oxygen is a good approximation, provided that the values of D
s
and h
s
are as derived in section A of the Appendix (additional file 1). Nevertheless, for completeness in the simulations of chemotherapy that follow we simulate the response of the full system (11)-(15).
Full model - the effect of chemotherapy
Given that the maximum rate of occlusion, δ, determines the intrinsic dynamics of our system (once the other parameters have been fixed), below we investigate how vessel status (δ low or high) influences therapeutic response. We simulate the response of different tumours to chemotherapy, comparing cases for which the system dynamics behind the invading front prior to treatment result in (i) a stable equilibrium, (ii) regular waves, and (iii) irregular waves. The only parameter change that is varied between these three cases is δ.
In Figure 6 we present typical snapshots showing how the tumour cell density varies across the spatial domain when chemotherapy is first applied (at t = 200, in dimensionless units, in all three cases). From panel A we estimate the speed of tumour invasion to be 0.002 cm/day; this value is of the same order of magnitude as the speed of 0.009 cm/day that has been observed in gliomas [30]. We note that the spatial extents, and thus the speeds, are similar in all three cases. Linearising ahead of the wave, the p-equation decouples and a travelling wave analysis indicates that there is a minimum wave speed which is independent of δ. Calculated speeds are in close agreement with the predicted minimum speeds.
A way of measuring the therapy's effect on tumour size is to monitor the tumour's spatial extent. Here we investigate how far the invading wave front has spread by determining the length of the interval in which the tumour cell density exceeds a threshold value, p
crit
, i.e. determining x
max
so that ∀x ∈ [x
max
, 1]: p(x, t) <p
crit
. Various imaging modalities, e.g. CT and MRI, function in this way; the tumour can only be detected wherever its density is above the "threshold of detection" (thus the real spatial extent of the tumour is in fact larger than that revealed by the medical image) [30]. In all simulations of chemotherapy we observed that therapy induces a similar response, with an initial decrease in x
max
, meaning that the therapy causes the invading wave front to recede.
In our simulations we observed that the therapy had a more pronounced effect on tumour burden than on the tumour's spatial extent. Therefore in what follows we only monitor the tumour burden, p
tot
, where p
tot
is defined by the Trapezoid rule
where p
i
= p(x
i
, t), x
i
= iΔx for i = 0, ..., n and Δx = L/n. In Figure 7 we show how the response to both a single bolus of chemotherapy and multiple boluses depends on the underlying system dynamics. We note that the initial tumour burden, p
tot
, is highest in the case of stable tumour dynamics (Figure 7, panel A) and lowest in the case of irregular oscillations (Figure 7, panel C), whereas the initial spatial extents are similar in all three cases (see Figure 6). Since occlusion is more pronounced (δ being high) when the dynamics are oscillatory, the varying initial values of p
tot
reveal that stable/strong vessels promote tumour growth more effectively than unstable/weak ones. Experimental observations of ovarian carcinoma spheroids have yielded similar trends: fast-growing tumours, that increased exponentially in volume, exhibited only small variations in the density of their vasculature, which mainly consisted of mature vessels [6]. Slow-growing tumours, on the other hand, had a larger proportion of immature vessels, and exhibited larger fluctuations in both their growth rate and their vessel density [6]. In Figure 7 we observe that in all cases initially therapy decreases the tumour burden. By comparing the tumour regrowth that occurs as the drug degrades, striking differences between the different cases are revealed. When the underlying dynamics of the control are such that the system evolves to a stable equilibrium (Figure 7, panel A), the cell population increases monotonically after a single bolus of chemotherapy. In this case the effect of therapy is unambiguous: the treated tumour is always smaller than the untreated one. Furthermore, multiple injections (corresponding to one bolus every 14 days in dimensional parameters and four times as much drug being administered compared to a single dosage) are clearly more effective at controlling tumour regrowth. When the untreated tumour undergoes regular or irregular oscillations (Figure 7, panels B and C, respectively), the oscillations during the regrowth phase are of larger amplitude than for the drug-free controls, making the effect of therapy harder to evaluate and less predictable. For example, in the case of irregular oscillations multiple injections do not provide an obvious improvement compared to only one round of therapy. The increase in the amplitudes of the oscillations that the therapy induces stems from the inherent oscillatory dynamics. Specifically, the initial decrease in the tumour cell density allows the vessels to recover from occlusion throughout the spatial domain. As the therapy wanes any tumour cells that remain have an abundant supply of oxygen. Therefore, the recovering tumour cell population attains higher total cell numbers than it does when no therapy is applied. The large cell population induces more extensive vessel occlusion, causing the tumour cell population to decline again to low cell numbers. In this way the large amplitude oscillations are sustained. We remark that in panel B, where the oscillations are regular, the period of the oscillations corresponds to roughly 17 days in dimensional terms. Slower/faster oscillations can be achieved by increasing/decreasing the maximum rate of vessel proliferation, η0, and the maximum rate of vessel occlusion, δ (ensuring δ/η0 > 1). Simulations with slower/faster oscillations (and others with D
v
≠ 1) yielded similar results as those presented above.