Oscillatory dynamics in a model of vascular tumour growth - implications for chemotherapy
© Stamper et al; licensee BioMed Central Ltd. 2010
Received: 12 February 2010
Accepted: 20 April 2010
Published: 20 April 2010
Investigations of solid tumours suggest that vessel occlusion may occur when increased pressure from the tumour mass is exerted on the vessel walls. Since immature vessels are frequently found in tumours and may be particularly sensitive, such occlusion may impair tumour blood flow and have a negative impact on therapeutic outcome. In order to study the effects that occlusion may have on tumour growth patterns and therapeutic response, in this paper we develop and investigate a continuum model of vascular tumour growth.
By analysing a spatially uniform submodel, we identify regions of parameter space in which the combination of tumour cell proliferation and vessel occlusion give rise to sustained temporal oscillations in the tumour cell population and in the vessel density. Alternatively, if the vessels are assumed to be less prone to collapse, stable steady state solutions are observed. When spatial effects are considered, the pattern of tumour invasion depends on the dynamics of the spatially uniform submodel. If the submodel predicts a stable steady state, then steady travelling waves are observed in the full model, and the system evolves to the same stable steady state behind the invading front. When the submodel yields oscillatory behaviour, the full model produces periodic travelling waves. The stability of the waves (which can be predicted by approximating the system as one of λ-ω type) dictates whether the waves develop into regular or irregular spatio-temporal oscillations. Simulations of chemotherapy reveal that treatment outcome depends crucially on the underlying tumour growth dynamics. In particular, if the dynamics are oscillatory, then therapeutic efficacy is difficult to assess since the fluctuations in the size of the tumour cell population are enhanced, compared to untreated controls.
We have developed a mathematical model of vascular tumour growth formulated as a system of partial differential equations (PDEs). Employing a combination of numerical and analytical techniques, we demonstrate how the spatio-temporal dynamics of the untreated tumour may influence its response to chemotherapy.
This manuscript was reviewed by Professor Zvia Agur and Professor Marek Kimmel.
Experimental observations suggest that blood vessels in solid tumours may be affected by the mechanical stress exerted on them by tumour cells. Immature vessels lacking pericytes are especially likely to be structurally unstable . As a tumour grows, the burden from tumour cells may cause the diameters of immature vessels in the tumour's interior to decrease [1, 2]. Indeed, completely collapsed, or occluded, vessels - those with a closed lumen - are found more frequently at higher tumour cell densities [1, 2]. On the other hand, if tumour cells are killed with cytotoxic therapy, the pressure on the vessels reduces, and vessel recovery may occur: the diameter of the vessels increases and the number of occluded vessels decreases [1, 2].
The model we develop in this paper is of reaction-diffusion type and describes a vascularised tumour growing in a one-dimensional spatial domain. We use the model to investigate how tumour cells and blood vessels may interact and influence the overall tumour growth dynamics. We assume that the tumour relies on oxygen delivered via the vasculature for sustained growth. The vascular density increases due to tumour angiogenesis, while it decreases due to vessel destabilisation and occlusion at a rate which is an increasing function of the ratio of tumour cell density to vessel surface area (reflecting the experimental observations mentioned above). Previously, other authors have used mathematical models to study vessel occlusion in tumours. As in our model, in the multiphase model developed by Breward et al.  angiogenesis and vessel occlusion act as opposing forces that regulate vessel growth, with vessel occlusion being switched on smoothly when the pressure exerted on the vessels by the cells exceeds a critical value. Although their model couples vessel and tumour cell dynamics, no oscillatory behaviour was reported - only stable growth dynamics were observed . An alternative approach to modelling vessel occlusion was proposed by Araujo and McElwain, who developed a continuum-mechanical model in which vascular collapse is linked to the evolution of residual stress within the growing tumour . In this model interactions between the tumour and (the implicit) vasculature (which is represented by a prescribed nutrient distribution) may generate oscillations in the tumour radius, albeit of very small amplitude .
Analysis of the model we derive in this paper reveals that the interplay between a tumour and its blood vessels may lead to strongly oscillatory dynamics. In particular, oscillations in the model variables arise if the vessels are structurally unstable, and thus are prone to occlusion (as opposed to stiffer and less compliant vessels which yield stable steady state solutions). The oscillations in our model are sustained by repeated cycles of tumour expansion and vessel occlusion, followed by episodes of tumour cell death and vessel recovery. As the tumour cells spread and colonise the spatial domain, regular or irregular spatio-temporal oscillations develop. These complex spatio-temporal dynamics may contribute to the reported spatio-temporal heterogeneity that characterises many vascularised solid tumours. Previously, oscillatory dynamics have been reported, for example by Arakelyan et al., who developed a spatially-averaged ordinary differential equation (ODE) model of tumour growth . Validation of the model was obtained from MRI data describing the growth of implanted human ovarian carcinoma . Furthermore, spatio-temporal irregularities have been observed in models which describe various aspects of the immune response to tumours [7–9].
Although our model does not explicitly account for blood flow, the highly irregular oscillations in vascular density and oxygen levels that the model generates are consistent with oxygen concentrations in tumours whose blood flow is chaotic. Chaotic blood flow is believed to be one of the main barriers to drug delivery , and the decreased vascular surface area and impaired blood flow of occluded vessels may also be detrimental to delivery [2, 11]. So as to gain insight into therapeutic implications, we investigate the impact of chemotherapy on the tumour. We focus on the cytotoxic drug doxorubicin, a widely used anti-cancer agent which may be used in free form or encapsulated in liposomes [12, 13]. Several models of doxorubicin administration have been developed [12, 14, 15]. For example, in  a compartmental model is used to compare the efficacy of free doxorubicin delivery versus liposome-bound delivery. The model in  distinguishes between vascular (plasma), extracellular and intracellular doxorubicin concentrations, accounting for the fact that the drug must enter the cells to be effective. A similar approach is used in  where cell kill depends on the amount of doxorubicin trapped within the cells. However, neither model accounts for both tumour cell proliferation and vascular remodelling, thus neglecting any effects that interactions between these phenomena are likely to induce.
On the contrary, our investigation concentrates on applying chemotherapy in different tumour growth scenarios, and shows how the intrinsic system dynamics (oscillatory or non-oscillatory, corresponding to compliant or non-compliant vessels, respectively) may alter treatment outcome. Specifically the ways in which the system recovers from therapy are found to vary markedly. This is consistent with simulations in  where the dynamics of tumour growth were modelled using a multiscale hybrid cellular automaton, and where the system's therapeutic response was influenced by the vessels' maturation status.
We organise this paper as follows. In Methods we develop our model of vascular tumour growth. In Results we show how, in a spatially homogeneous version of our model (in which spatial effects are neglected and fast kinetics are assumed for oxygen), the tumour cell population and the vessel density evolve to a stable co-existence equilibrium, or, if the steady state is unstable, undergo oscillations (the latter occurring when the maximum rate of vessel occlusion exceeds a critical value). Then, by reintroducing spatial effects, we investigate the way in which a small population of tumour cells seeded at one end of the domain invades the tissue. In agreement with analysis of other reaction-diffusion systems, we find that the form of the solutions behind the invading tumour cell front depends on the stability of the spatially homogeneous co-existence steady state behind the front . Specifically, when the co-existence steady state is stable, the system evolves to this stable equilibrium in the wake of the invading front. In contrast, in the case of an unstable co-existence steady state, the system is subject to either regular or irregular oscillations. We are able to predict the type of oscillations that will arise for particular model parameters by analysing a system of λ-ω type that approximates the dynamics close to the unstable steady state; this approach has previously been applied by Sherratt . We then investigate the effect of chemotherapy. We conclude with a discussion of our results, and suggestions for future research.
where J u denotes the flux and f u represents the net sources and sinks of species u. In what follows we adapt equation (1) for our dependent variables.
Tumour cells, p(x, t)
where β is the maximum proliferation rate per unit of cell density and s β is the level of oxygen at which the proliferation rate is half-maximal. We assume that cell death occurs naturally due to apoptosis at rate d p p, where d p is a positive constant. To account for tumour cell death induced by the cytotoxic drug, we assume that the drug kills proliferating cells in a concentration-dependent manner  and that the maximum rate of cell kill is equal to the rate of proliferation (a similar assumption for the action of a cytotoxic drug was used in ). We note that c = 0 implies that the rate of cell proliferation is unaffected (i.e. all cells are proliferating and there is no drug-induced cell kill), c = K c implies zero net contribution from proliferation (i.e. proliferation has ceased in half the cell population, which instead dies), and c → ∞ implies no contribution from proliferation (i.e proliferation has ceased in all cells, which instead die). The precise mechanism by which doxorubicin is anti-proliferative is not known, but it is in part due to intercalation into DNA .
Regarding tumour cell movement, we allow for random motion of the cells and so specify , where D p is the random motility co-efficient which is assumed to be constant.
where K c is a positive constant denoting the concentration at which the cell kill is half-maximal and D p , β, s β and d p are positive constants.
Tumour vasculature, v(x, t)
While not explicitly introducing VEGF or any other angiogenic factor into our model, we assume that the tumour stimulates angiogenesis in such a way that the vasculature undergoes logistic growth characterised by a constant rate of proliferation, η0, and a carrying capacity, V0. We note here that a more realistic model of tumour-induced angiogenesis could include VEGF (e.g. VEGF could be secreted by the tumour cells, especially when the level of oxygen is low). The level of angiogenesis could then be linked explicitly to the concentration of VEGF, for example by allowing η0 to be an increasing, saturating function of VEGF. However, since our goal is to gain insight into the dynamics that arise from tumour cell proliferation and vessel occlusion, we refrain from adding this extra complexity and assume that η0 is constant (equivalent to VEGF being maintained at a constant level).
where the constant K δ denotes the value of p/v at which the rate of occlusion attains its half-maximal value, δ/2.
To account for movement of the endothelial cells (ECs) that constitute the vessels we include random motion, assuming that this occurs with a constant diffusivity co-efficient, D v . Since VEGF is not explicitly included in the model we neglect EC movement due to chemotaxis.
where D v , δ, K δ , η0 and V0 are positive constants.
Oxygen level, s(x, t)
We assume that the concentration of oxygen, s(x, t), increases due to delivery from the vasculature. Transcapillary exchange of molecules from the blood occurs mainly by diffusion and convection . For simplicity we suppose that diffusive effects dominate and that the diffusive flux across the vessel wall is proportional to the vascular surface area and to the difference between the concentrations of the substance in the plasma and in the interstitium . Thus we model the delivery of oxygen from the blood as h s (s b - s)v, where h s is the vascular permeability constant , s b is the (assumed constant) concentration of oxygen in the blood, s is the concentration of oxygen in the tissue and v is the vascular surface area per unit tissue volume as explained above. We note that as the oxygen flows through the tissue the level of oxygen in the blood, s b , decreases along each vessel due to the diffusive flux into the tissue. However, since we do not account for vessel morphology we make the simplifying assumption that s b is a constant.
After transcapillary transfer, the transport of oxygen within the interstitium occurs by diffusion and convection . Neglecting the convective part for simplicity, we assume that the transport of oxygen within the tissue is via diffusion only so that .
We assume that oxygen is consumed by tumour cells, at rate σ p ps, and by other cellular species (e.g. fibroblasts and macrophages) resident in the tissue, at rate d s s. An alternative form for tumour oxygen consumption could be obtained by assuming that the consumption rate is of Michaelis-Menten type. By assuming that d s > 0 we ensure that, in the absence of tumour cells, the level of oxygen evolves to a finite, tissue-specific equilibrium value. The particular value of the parameter d s , and the associated equilibrium oxygen value, will be tissue-specific and could vary as the tumour and the surrounding stroma evolve.
where D s , h s , s b , σ p and d s are positive constants.
Plasma level of drug, c b (t)
so that the drug decays exponentially after being injected. In the case of multiple injections we assume that the doses are such that the plasma concentration regains its initial value at each injection (achieved mathematically in (6) by making the time variable modulo T per ). More precisely, we assume T1/2 ≪ T per , as is the case for doxorubicin administration, where the half-life is roughly one day and the period is on the order of weeks .
Concentration of drug in tumour tissue, c(x, t)
where D c , h c and d c are positive constants. In equation (7) the vascular delivery term is similar to that introduced for oxygen in (4), while the decay term represents drug uptake by both tumour cells and other cells present in the extracellular matrix.
Boundary and initial conditions
where p i is an arbitrary constant and s i = h s s b v/(h s v + σ p p + d s ). This implies that the tumour cells are introduced at the left-hand boundary of a tissue in which the vasculature is at its carrying capacity, V0, and the oxygen level is at quasi-steady state. In addition there is no drug in the tissue initially. Before analysing the model it is convenient to nondimensionalise it as explained in the next section.
are dimensionless parameters. A discussion of the parameter values that we use is contained in section A of the Appendix (additional file 1).
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 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.
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
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 . 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 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 . 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) . 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  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 . 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 . 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].
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 δ.
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) . 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 this paper we have developed a PDE model which describes tumour growth within a vascularised tissue where the tumour cells rely on the vasculature for their supply of oxygen. Analysis of a spatially homogeneous submodel, for which the oxygen profile is assumed to be quasi-steady, showed that the tumour cell and vessel densities either evolve to a stable equilibrium or undergo oscillations. We found that the latter is more likely to occur if the vessels are prone to occlusion (δ being large, corresponding to immature vessels) and that the oscillations arise from interactions between the tumour cells and the compliant blood vessels.
Oscillations in mathematical models of tumour growth have been described by others [5, 31]. However, different mechanisms are responsible for those oscillatory dynamics. In  the oscillations are due to time-delays in an ODE-system, while in  the underlying cause is vascular adaptation induced by VEGF secretion by quiescent tumour cells. To our knowledge, this is the first model of vascularised tumour growth in which the direct interactions between the tumour cells and the vasculature lead to sustained oscillations of significant amplitude (very small oscillations were observed in the model developed by Araujo and McElwain ). Experimental investigations indicate that tumour cells have the capacity to occlude vessels through the pressure they exert on the vessel wall [1, 2]. Furthermore, the observation that vessel occlusion is reversible (i.e. if the tumour burden is relieved by killing off tumour cells cytotoxically, then the vessel diameter increases [1, 2]) leads us to believe that oscillations of the kind we report here are physically realistic and may thus contribute to spatio-temporal heterogeneities in solid tumours. At present we lack clinical data to determine whether such oscillations are clinically observable. However, our simulations (see Figure 7) suggest that oscillations could indeed be detected from spatially-averaged clinical data, provided that the data is adequately resolved in time.
Our analysis of tumour invasion due to random cell movement revealed varying spatio-temporal patterns depending on the stability of the spatially homogeneous co-existence steady state behind the invading front. Following invasion in the case of a stable equilibrium, the system settled to this steady state. In the case of an unstable equilibrium and oscillatory solutions in the spatially homogenous system, spatial effects induced either regular or irregular spatio-temporal oscillations. By approximating our system by one of λ-ω form near the Hopf bifurcation, we were able to predict the type of oscillation that will occur. This technique has been applied previously to predator-prey models where similar spatio-temporal oscillations are observed .
When investigating tumour response to chemotherapy we focussed on the different growth scenarios that can arise as a result of tumour invasion. We simulated chemotherapy in the form of either a single bolus or multiple, consecutive boluses; other possible modes of delivery include continuous infusion, investigated e.g. in . Our numerical simulations reveal that following administration of a single bolus of doxorubicin to a tumour which, in the absence of therapy, is growing toward a steady state behind the invading front, a reduction in tumour burden is achieved. During the regrowth period that occurs while the drug is degrading, the tumour increases monotonically in size, but the size of the tumour cell population is smaller than that of the untreated control. If the tumour undergoes regular or irregular oscillations behind the invading front, the initial response is also a decrease in the tumour burden. However, during the regrowth phase, the size of the tumour cell population fluctuates strongly, more so than without therapy, and the tumour cell population attains both higher and lower values than for the drug-free cases. The therapy reinforces the underlying cause of the oscillatory dynamics (the relief of occlusion due to pressure from the expanding tumour), leading to an increase in the amplitude of the oscillations. Not only would such irregular spatio-temporal dynamics make the therapeutic effect harder to evaluate (since the timing of measurements would be crucial), but they may also explain why chemotherapy sometimes fails. In particular, our model suggests that clinical data showing short-time regression and longer-time relapse may be the result of tumours whose dynamics are oscillatory. However, in order to generate quantitative predictions, comprehensive and systematic experiments in which all model parameters can be estimated from the same cell line are urgently needed. In the absence of such data, a detailed sensitivity analysis (where system parameters are varied within physiological ranges) could provide additional information about the relative importance of individual parameters on the system dynamics.
It is interesting to compare our simulations of a single bolus of doxorubicin to the simulations presented in , wherein a modified version of the multiscale model developed in  was studied. The modifications in  relate specifically to the effect that tumour cells have on the pre-existing vasculature. Vessel remodelling is incorporated at each time step so that vessels surrounded by tumour cells are considered immature/destabilised, whereas those without tumour cells are considered mature. The inclusion of vessel destabilisation/maturation has pronounced effects on the therapeutic response of the system in : after a single bolus, in the case where the vessels undergo destabilisation/maturation, pronounced oscillations in the tumour cell population occur as the tumour cells recover from the therapy; when destabilisation/maturation is neglected, the tumour cell recovery is more regular and the tumour cell population increases steadily. These dynamic responses are similar to the results we obtain when the tumour is exposed to a single bolus (Figure 7): during the regrowth phase, the tumour undergoes oscillations (when the underlying equilibrium is unstable) or continuously increases (when the underlying equilibrium is stable). We remark here that in our model the change from monotonically increasing tumour recovery to oscillatory tumour recovery occurs as the maximum strength of vessel occlusion is increased. Therefore, our oscillatory vessels may be regarded as more unstable than those present in tumours that recover without oscillations, a phenomenon that would be consistent with the ideas presented in  that vessel destabilisation introduces an irregular tumour growth pattern due to blood flow heterogeneities.
As with any mathematical description of a biological system, modifications could be made to increase the realism of our model. For example, we could investigate the effects that the angiopoietins, Ang-1 and Ang-2, have on the developing vasculature. To do this we would need to separate the vessels into two subpopulations consisting of immature vessels (destabilised by Ang-2) and mature vessels (stabilised by Ang-1). Immature vessels may be particularly sensitive to occlusion since they typically lack a pericyte covering . Another obvious extension involves introducing VEGF explicitly into the model, allowing the rate of EC proliferation to be a saturating function of the VEGF concentration. We anticipate that these model extensions would yield similar qualitative behaviour, e.g. if the maximum rate of vessel occlusion (δ) exceeded the rate of EC proliferation (η0) then the oscillatory solutions would persist. By introducing VEGF it would also be possible to include chemotaxis with vessels moving up spatial gradients of VEGF. The inclusion of chemotaxis is likely to affect the spatio-temporal patterns that emerge as a result of tumour invasion.
Regarding the way in which we have modelled chemotherapy it is also possible to amend the model to more accurately describe doxorubicin activity. For example, we could investigate the effects of changing the functional form used to model cell kill in equation (2). As found in , the model predictions may depend on the numerical value of the threshold at which the rate of cell kill switches from low to high (in our model, this value is represented by the parameter K c ). We could also investigate the tumour's response to therapy if the drug targets both tumour cells and proliferative ECs.
Despite the simplicity of the model, the simulations of doxorubicin administration presented in this paper still give considerable insight into how the response to chemotherapy may be strongly influenced by the underlying system growth dynamics. A possible direction for further research involves extending our model to distinguish between proliferating and quiescent tumour cells in order to determine how the dynamics are affected by cell heterogeneity. We could also extend our model to investigate the impact that vessel normalisation strategies, designed to increase the stability of the immature tumour vessels, have on systems with irregular spatio-temporal dynamics. Possible candidates for vascular-targetting agents include the anti-angiogenic drug DC101 which transiently increases pericyte recruitment, hence increasing vessel stability, oxygenation, and drug delivery . Modelling vascular-targetting treatments could be achieved by making the occlusion parameter δ (which appears in equation (3)) a decreasing function of the vascular-targetting drug concentration. We anticipate that this could render the system more stable, possibly enhancing the therapy's performance.
In this paper we have developed a deterministic model of vascular tumour growth formulated as a system of reaction-diffusion equations. Using a combination of numerical methods and analytical techniques we showed how the system's spatio-temporal dynamics (and, in particular, tumour-invasion patterns) are influenced by the dynamics of the associated spatially-uniform system. Depending on the interplay between the proliferating tumour cells and the collapsible vessels, we observed patterns of tumour invasion ranging from regular travelling waves to waves with highly irregular spatio-temporal structure. The latter were characterised by highly responsive, compliant vessels while the former were characterised by stiffer, less compliant vessels. Simulations of chemotherapy illustrated how treatment outcome depends on the dynamics of the untreated tumour. In particular we found that the therapeutic efficacy could be dramatically reduced by the presence of immature vessels which contribute to the emergence of oscillatory dynamics.
Reviewer's report 1
Prof. Zvia Agur, The Institute for Medical BioMathematics (IMBM), Israel
This paper focuses on vessel occlusion. Analysis of a mathematical model for vessel occlusion, in the form of a system of reaction-diffusion equations, shows that increased occlusion may give rise to sustained oscillations in tumor cell density and in vessel density. The paper suggests that oscillations in tumor size may compromise the efficacy of Doxorubicine chemotherapy. I would like to mention a few points which may be considered in order to improve the paper.
1. The basic model deliberately simplifies tumor growth and vascularisation dynamics, and couples these simple dynamics through vessel occlusion (see, e.g., Fig 2). However, occlusion is not a frequent phenomenon in tumor angiogenesis, where angiogenesis and tumor growth are always coupled through various other factors, nor is it prerequisite for oscillations in tumor and vascular dynamics (see below).
2. The paper makes some claims, which require substantiation, or should be reconciled with literature, or should be altogether eliminated. In the Discussion section (P. 15, Paragraph 3) the authors state: "To our knowledge, this is the first model of vascularised tumor growth in which the direct interactions between the tumor cells and the vasculature lead to sustained oscillations of significant amplitude". However, I believe this statement should be eliminated. Both some of the papers cited in the present manuscript and others have already showed the same result. For example, in , we have introduced and analysed several alternative general tumor angiogenesis models (with/without vascular maturation or with/without the effects of growth factors). Our results show that whenever a time-delay is introduced into the tumor proliferation or the neo-vascularization process, Hopf points are found, leading to oscillatory behavior. While it is recognized that time-delay will often elicit Hopf points, in , it was shown that the latter were to be found in any angiogenesis model with time-delay. There are other examples of angiogenesis models showing oscillatory tumor and vasculature dynamics as a result of coupling between these two processes, for example, as a result of certain ratios between growth factors' reaction rates. Suggestion: The authors should crystallize the novelty in their results and be precise about the conditions for obtaining it. In this context, I believe they should make a distinction between occlusion and hypoxia, which may be caused by many other factors.
3. Another example. On page 3 parag 2 Stamper et al write: " Analysis of the model we derive in this paper reveals that the interplay between tumor and its blood vessels may lead to strongly oscillating dynamics... if the vessels are structurally unstable". Previous models with no occlusion assumption have shown the same result. Some of them, but not all, are cited in the present manuscript (e.g., ).
4. One of the main conclusions of the paper is that ("Conclusions", p. 21, Parag 2): "Simulations of the chemotherapy illustrated how treatment outcome depends on the dynamics of the untreated tumor...". I believe this is already known to any scientist involved in cancer research. In any case it has been shown in many papers. I think this sentence should be suppressed.
5. I suggest that sentences such as (Page 18, parag 2) "We found that ...oscillations arise from the interactions between the tumor cells and the compliant blood vessels" would rather be eliminated.
6. I have not found any reference to the pharmakodynamics of doxorubicine. What is the assumed drug action mechanism? Why does the model represent doxorubicine and no other chemotherapy (half-life)?
Top sum up, I think that the paper will benefit from a careful editing where the main messages are spelled out consistently through the Abstract, Results, and Conclusions sections, and where the special assumptions and simplifications are recognised so that conclusions and implications are more restricted.
We agree with Professor Agur that our model is not the first one of vascular tumour growth to produce oscillatory dynamics - that is why we cite the seminal work of Agur and coworkers, amongst others. The essential novelty in our paper arises from the investigation of how the interplay between the oscillations associated with the underlying time-dependent model and diffusive effects can give rise to highly irregular spatio-temporal dynamics within the tumour mass and their impact on the tumour's response to therapy.
Reviewer's report 2
Prof. Marek Kimmel, Department of Statistics, Rice University, USA
A comment concerning the catalogue of spatial oscillatory patterns. In the paper  by Marciniak-Czochra and Kimmel, a new model of field carcinogenesis, is explored inspired by lungcancer precursor lesions, which includes dynamics of a spatiallydistributed population of pre-cancerous cells c(t; x), constantly supplied by an influx of mutated normal cells. Cell proliferation is controlled by growth factor molecules bound to cells, b(t; x). Free growth factor molecules g(t; x) are produced by precancerous cells and may diffuse before they become bound to other cells. The model shows an invading front of oscillatory solutions, which correspond to formation of multiple spatially isolated lesions of pre-cancerous cells and which become in the limit stable spike solutions. This effect is very different from the one obtained in Stamper et al. (although some analogies between oxygen and growth factor supply exist), among other because only one of the three equations features a diffusion term.
We thank Professor Kimmel for directing us to this interesting publication. There are some similarities with our work, in that we consider diffusion-driven instability (DDI) as a mechanism for destabilising spatially homogeneous solutions. However an important difference between our work and that of Marciniak-Czochra and Kimmel is that we also use lambda-omega analysis to investigate the stability of time-periodic model solutions to spatial perturbations. In this way we provide a possible explanation for the highly irregular spatio-temporal dynamics that characterise many vascularised tumours and show that such dynamics may have a deleterious impact on the tumour's response to chemotherapy.
IJS gratefully acknowledges the support of the European Commission under MCRTN. The views expressed here are not necessarily those of the European Commission. PKM was partially supported by a Royal Society-Wolfson Merit Award and NIH grant U56CA113004 from the National Cancer Institute. As a one-time exception to the publishing policy of Biology Direct the articles in this series are being published with two reviewers.
- Padera TP, Stoll BR, Tooredman JB, Capen D, di Tomaso E, Jain RK: Cancer cells compress intra-tumour vessels. Nature. 2004, 427: 695-10.1038/427695a.PubMedView Article
- Griffon-Etienne G, Boucher Y, Brekken C, Suit HD, Jain RK: Taxane-induced apoptosis decompresses blood vessels and lowers interstitial fluid pressure in solid tumours: clinical implications. Cancer Res. 1999, 59: 3776-3782.PubMed
- Breward CJW, Byrne HM, Lewis CE: A multiphase model describing vascular tumour growth. Bull Math Biol. 2003, 65: 609-640. 10.1016/S0092-8240(03)00027-2.PubMedView Article
- Araujo RP, McElwain DLS: The role of mechanical host-tumour interactions in the collapse of tumour blood vessels and tumour growth dynamics. J Theor Biol. 2006, 238: 817-827. 10.1016/j.jtbi.2005.06.033.PubMedView Article
- Arakelyan L, Vainstein V, Agur Z: A computer algorithm describing the process of vessel formation and maturation, and its use for predicting the effects of anti-angiogenic and anti-maturation therapy on vascular tumour growth. Angiogenesis. 2002, 5: 203-214. 10.1023/A:1023841921971.PubMedView Article
- Arakelyan L, Merbl Y, Agur Z: Vessel maturation effects on tumour growth: validation of a computer model in implanted human ovarian carcinoma spheroids. Eur J Canc. 2005, 41: 159-167. 10.1016/j.ejca.2004.09.012.View Article
- Matzavinos A, Chaplain MAJ: Travelling-wave analysis of a model of the immune response to cancer. C. R. Biologies. 2004, 327: 995-1008. 10.1016/j.crvi.2004.07.016.PubMedView Article
- Matzavinos A, Chaplain MAJ: Mathematical modelling of the spatio-temporal response of cytotoxic T-lymphocytes to a solid tumour. Math Med Biol. 2004, 21: 1-34. 10.1093/imammb/21.1.1.PubMedView Article
- Owen MR, Sherratt JA: Pattern formation and spatiotemporal irregularity in a model for macrophage-tumour interactions. J Theor Biol. 1997, 189: 63-80. 10.1006/jtbi.1997.0494.PubMedView Article
- Jain R: Determinants of tumor blood flow: a review. Canc Res. 1988, 48: 2641-2658.
- Jain RK: Normalization of tumor vasculature: an emerging concept in antiangiogenic therapy. Science. 2005, 307: 58-62. 10.1126/science.1104819.PubMedView Article
- El-Kareh AW, Secomb TW: A mathematical model for comparison of bolus injection, continuous infusion, and liposomal delivery of doxorubicin to tumor cells. Neoplasia. 2000, 2: 325-338. 10.1038/sj.neo.7900096.PubMedPubMed CentralView Article
- Rahman A, Carmichael D, Harris M, Roh JK: Comparative pharmacokinetics of free doxorubicin and doxorubicin entrapped in cardiolipin liposomes. Canc Res. 1986, 46: 2295-229.
- Jackson TL: Intracellular accumulation and mechanism of action of doxorubicin in a spatio-temporal tumor model. J Ther Biol. 2003, 220: 201-213. 10.1006/jtbi.2003.3156.View Article
- Ribba B, Marron K, Agur Z, Alarcón T, Maini PK: A mathematical model of doxorubicin treatment efficacy for non-Hodgkin's lymphoma: investigation of the current protocol through theoretical modelling results. Bull Math Biol. 2005, 67: 79-99. 10.1016/j.bulm.2004.06.007.PubMedView Article
- Sherratt JA, Eagan BT, Lewis MA: Oscillations and chaos behind predator-prey invasion: mathematical artifact or ecological reality?. Phil Trans R Soc Lond B. 1997, 352: 21-38. 10.1098/rstb.1997.0003.View Article
- Sherratt JA: Periodic travelling waves in cyclic predator-prey systems. Ecology Letters. 2001, 4: 30-37. 10.1046/j.1461-0248.2001.00193.x.View Article
- Orme ME, Chaplain MAJ: A mathematical model of vascular tumour growth and invasion. Math Comput Modelling. 1996, 23: 43-60. 10.1016/0895-7177(96)00053-2.View Article
- Secomb T, Hsu R, Dewhirst MW: Synergistic effects of hyperoxic gas breathing and reduced oxygen consumption on tumor oxygenation: a thereoretical model. Int J Radiation Oncology Biol Phys. 2004, 59: 572-578.View Article
- Edelstein-Keshet L: Mathematical models in biology. 1988, NY: McGraw-Hill Inc
- Casciari JJ, Sotirchos SV, Sutherland RM: Variations in tumour cell growth rates and metabolism with oxygen concentration, glucose concentration and extracellular pH. J Cell Physiol. 1992, 151: 386-394. 10.1002/jcp.1041510220.PubMedView Article
- Fornari FA, Randolph JK, Yalowich JC, Ritke MK, Gewirtz DA: Interference by doxorubicin with DNA unwinding in MCF-7 breast tumor cells. Mol Pharm. 1994, 45: 649-656.
- Webb SD, Owen MR, Byrne HM, Murdoch C, Lewis CE: Macrophage-based anti-cancer therapy: modelling different modes of tumour targeting. Bull Math Biol. 2007, 69: 1747-1776. 10.1007/s11538-006-9189-2.PubMedView Article
- Jain RK: Transport of molecules, particles, and cells in solid tumors. Annu Rev Biomed Eng. 1999, 01: 241-263. 10.1146/annurev.bioeng.1.1.241.View Article
- NAG Fortran library routine document D03PCF/D03PCA. [http://www.nag.co.uk/numeric/Fl/manual/pdf/D03/d03pcf.pdf]
- Murray JD: Mathematical Biology. 1993, Berlin: Springer VerlagView Article
- Sherratt JA: Invading wave fronts and their oscillatory wakes are linked by a modulated travelling phase resetting wave. Physica D. 1998, 117: 145-166. 10.1016/S0167-2789(97)00317-5.View Article
- Guckenheimer J, Holmes P: Nonlinear oscillations, dynamical systems, and bifurcations of vector fields. 1983, NY: SpringerVerlagView Article
- Malchow H, Petrovskii SV: Dynamical stabilization of an unstable equiibrium in chemical and biological systems. Math Comp Mod. 2002, 36: 307-319. 10.1016/S0895-7177(02)00127-9.View Article
- Swanson KR, Bridge C, Murray JD, Jr ECA: Virtual and real brain tumors: using mathematical modeling to quantify glioma growth and invasion. J Neurological Sci. 2003, 216: 1-10. 10.1016/j.jns.2003.06.001.View Article
- Byrne HM, Owen MR, Alarcón T, Murphy J, Maini PK: Modelling the response of vascular tumours to chemotherapy: a multiscale approach. Math Mod Meth Appl Sci. 2006, 16: 10.1142/S0218202506001522.
- Lin MI, Sessa WC: Antiangiogenic therapy: creating a unique "window" of opportunity. Cancer Cell. 2004, 6: 529-531.PubMed
- Agur Z, Arakelyan L, Daugulis P, Ginosar Y: Hopf point analysis for angiogenesis models. Discrete and Continuous Dynamical Systems Series B. 2004, 4: 29-38. 10.3934/dcdsb.2004.4.29.View Article
- Arakelyan L, Merbl Y, Daugulis P, Ginosar Y, Vainstein V, Selitser V, Kogan Y, Harpak H, Agur Z: Multi-scale analysis of angiogenic dynamics and therapy. Cancer Modelling and Simulation. Edited by: Preziosi L. 2003, Chapman and Hall/CRC, 185-220.
- Marciniak-Czochra A, Kimmel M: Reaction-difusion model of early carcinogenesis: The effects of influx of mutated cells. Math Model Nat Phenom. 2008, 3: 90-114. 10.1051/mmnp:2008043.View Article
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.