Capturing the dynamics of a hybrid multiscale cancer model with a continuum model

Cancer is a complex disease involving processes at spatial scales from subcellular, like cell signalling, to tissue scale, such as vascular network formation. A number of multiscale models have been developed to study the dynamics that emerge from the coupling between the intracellular, cellular and tissue scales. Here, we develop a continuum partial differential equation model to capture the dynamics of a particular multiscale model (a hybrid cellular automaton with discrete cells, diffusible factors and an explicit vascular network). The purpose is to test under which circumstances such a continuum model gives equivalent predictions to the original multiscale model, in the knowledge that the system details are known, and differences in model results can be explained in terms of model features (rather than unknown experimental confounding factors). The continuum model qualitatively replicates the dynamics from the multiscale model, with certain discrepancies observed owing to the differences in the modelling of certain processes. The continuum model admits travelling wave solutions for normal tissue growth and tumour invasion, with similar behaviour observed in the multiscale model. However, the continuum model enables us to analyse the spatially homogeneous steady states of the system, and hence to analyse these waves in more detail. We show that the tumour microenvironmental effects from the multiscale model mean that tumour invasion exhibits a so-called pushed wave when the carrying capacity for tumour cell proliferation is less than the total cell density at the tumour wave front. These pushed waves of tumour invasion Tanvi V. Joshi Stuttgart Research Center Systems Biology, University of Stuttgart, 70569, Germany Tel.: +49 176 76 312458 E-mail: tanvij03@gmail.com Daniele Avitabile Centre for Mathematical Medicine and Biology, School of Mathematical Sciences, University of Nottingham, NG7 2RD, UK Markus R. Owen Centre for Mathematical Medicine and Biology, School of Mathematical Sciences, University of Nottingham, NG7 2RD, UK.

propagate by triggering apoptosis of normal cells at the wave front. Otherwise, numerical evidence suggests that the wave speed can be predicted from linear analysis about the normal tissue steady state.
Keywords Multiscale · Continuum · Wave speed · Pushed waves 1 Introduction Cancer is a complex disease not restricted to a single biological scale but involving processes occurring over multiple spatial scales, ranging from the subcellular scale (for example, progression through the cell cycle) to the tissue scale (for example, angiogenesis, vascular remodelling), and temporal scales (for example, oxygen diffusion on a relevant tissue length-scale occurs in minutes while the doubling time for cells is usually in days/weeks).
Many models developed previously focused on a single aspect of tumour growth, for example, cell cycle (Alarcón et al, 2004), cell apoptosis and necrosis (Byrne and Chaplain, 1998), diffusion of angiogenic factors (Chaplain and Stuart, 1991) or tumour angiogenesis (Balding and McElwain, 1985). Although these models provide a good insight into the underlying mechanisms of these processes, they fail to address issues of how different phenomena occurring at different scales couple together to influence tumour growth. Many have already developed models to address this issue of bridging scales. To name a few, (Ribba et al, 2006) presented a multiscale model including key genes, cellular and tissue level dynamics and radiosensitivity dependence on cell cycle phase to study tumour response to irradiation protocols, (Anderson et al, 2007) shows how three multiscale individual cell based models emerging from cancer related processes at extracellular, cellular and subcellular scales are related to one another and can be used to bridge the three spatial scales, (Stolarska et al, 2009) developed multiscale models incorporating the mechanical processes involved in cell and tissue movement and studied them in the context of tumour growth, (Alarcón et al, 2005) developed a multiple scale model accounting for interactions between normal and cancerous cells, cell division and apoptosis, transport of oxygen and VEGF, blood flow and structural adaptation of vessels and coupling between these processes. The multiple scale model from (Alarcón et al, 2005) was later extended to include several features like cell movement (Betteridge et al, 2006), angiogenesis and vascular remodelling (Owen et al, 2009) and macrophage based gene therapy .
In this paper, we develop a model for vascular tumour growth, based on the multiscale model (MM) from , but consisting of only Partial Differential Equations (PDEs). The purpose of developing this continuum model (CM) is to test under which circumstances will such a model will be equivalent to the multiscale model and to highlight the advantages and the disadvantages of the two types of model. Comparing models based on the same set of features is beneficial since in such cases all the system details are known, and, any differences in the simulation results can be explained in terms of the model features. In contrast, differences in the simulation results of a model and the corresponding experimental data can be difficult to explain due to presence of several unknown factors. The process of developing a continuum model based on the dynamics observed in the multiscale model uncovers some important but not apparent features of the multiscale model. Additionally, continuum models enable us to perform certain analyses, like determining the steady states of the system and investigating the effect of model parameters and initial conditions on the evolution of the model variables.

Model development
The model we develop is a continuum model based on the multiscale model by Owen et al (2011). The multiscale model is based on a hybrid cellular automaton concept. The model comprises three layers, the cellular layer, the diffusibles like oxygen and growth factors (VEGF) and the vasculature describing the formation and function of a vessel network. We illustrate the three layers, the coupling between them and the order of computation of all the processes, occurring at various length scales, over a time step in the simulation in the flowchart in Fig. 1. Note that, in this article, we only study normal tissue expansion and tumour invasion. Hence, we do not include macrophages and processes associated with it from Owen et al (2011).  Owen et al (2011) To account for the three layers of the multiscale model and the processes associated with them, we define the following variables in the continuum model, as functions of the vector of space co-ordinates x and time t, -Oxygen Concentration: -VEGF Concentration: -Normal Cell Density: -Tumour Cell Density: -Immature Tip Cell Density: -Mature Tip Cell Density -Blood Vessel Surface Density: C(x,t) mmHg V (x,t) nM n(x,t) cells/mm 3 r(x,t) cells/mm 3 e(x,t) cells/mm 3 m(x,t) cells/mm 3 b(x,t) mm 2 /mm 3 In the remainder of this section, we will describe in detail the implementation of the three layers in the multiscale model and the corresponding equations in the continuum model.

Cell population
The cellular layer in the multiscale model accounts for individual cells of different types and certain subcellular processes associated with them. The layer is modelled as a lattice in either two or three spatial dimensions. However, here we focus only on two spatial dimensions. Each lattice site can hold a variable number of cells, depending on the cell type. At each time step, the cells progress through their respective cell cycles and divide into daughter cells if appropriate conditions are met, move within the lattice and/or die depending on the external conditions, like nutrient availability.
In the continuum model we consider cell densities, which correspond to the number of cells per unit volume in the multiscale model. The cell density changes according to the rate of movement, proliferation and death of the cells.
Below, we discuss in detail the implementation of cell movement, proliferation and death in the multiscale model and derive a general form for each of these processes and for each cell type. The PDE for a given cell type will then be derived combining all the corresponding relevant terms.

Cell movement
In the multiscale model, cells perform a random walk within the lattice. Some cells, like the tip cells sprouting from pre-existing blood vessels, move randomly, biased by the VEGF gradient. This type of movement is termed as chemotaxis.
In 2D, the probability, P(x, y,t) of a cell moving from site x to site y, a site in the Moore neighbourhood of x (i.e. the eight connected neighbours, including the diagonal sites), in time ∆t is given as where N(x,t) is the total number of cells at site x, V (x,t) is the VEGF concentration at site x, D cell is the random motility coefficient for cell type, N cell m is the carrying capacity, χ cell is the chemotactic sensitivity and d x,y is the distance between sites x and y. The probability of a cell not moving is one minus the sum of probabilities of moving to the eight sites in the Moore neighbourhood. Eq. (1) is used for all cell types and the superscript cell denotes that the parameter value depends on the cell type, like normal, tumour or tip cells, attempting to move. Also, random movement without chemotaxis is obtained by setting χ cell = 0 (e.g. for normal and tumour cells). Note that, when N(y,t) > N cell m or VEGF gradients are sufficiently steep, Eq. (1) can give negative values. In the former case there is no space available at site y, and in the latter case chemotaxis inhibits movement to that site, so we set the probability P(x, y,t) to 0 in such cases.
When χ cell = 0, a simple constraint on the time step, ∆t ≤ ∆ x 2 /(3D cell ) (where ∆ x is the lattice spacing), ensures that the sum of movement probabilities (and hence all individual probabilities) cannot exceed one. However, when χ cell = 0, steep VEGF gradients can violate this. Thus, whenever the sum of the probabilities exceeds one, we divide each probability by that sum. Correspondingly, the probability of not moving becomes zero. We note that this requirement arises rarely in the multiscale model simulations, and this normalisation was not previously implemented in e.g. Owen et al (2011).
Eq. (1) is similar to the volume filling model developed by Painter and Hillen (2002) and its continuum limit, for cell density, u(x,t), is given as (2) Note that the linear form of the volume filling term means that it has no effect on the diffusion term. Also, currently in the model, χ cell is non zero only for tip cells which have the highest carrying capacity. Hence u ≤ N cell m .

Cell proliferation
In the multiscale model, the phase φ of the cell cycle in each cell is modelled as an ordinary differential equation (ODE) where C cell φ is the oxygen concentration at which the speed is half maximal and T cell min is the minimum period of cell cycle. Eq. (3) is used for normal and tumour cell proliferation. φ is a number between 0 and 1, where 0 marks the beginning of the cell cycle and 1 marks the end of cell cycle after which the cell divides. Hence, when φ = 1, the cell divides and φ is (re)set to zero for the parent and daughter cells.
If C is constant, a new cell is born every T cell min (C cell φ +C)/C min, as this is the time taken by φ to go from 0 to 1. Thus, the cell density, u, increases as where u 0 is the cell density at t = 0.
Additionally, in the multiscale model, every cell type has a carrying capacity. This means that the cells do not divide into the target site if the total number of cells there is greater than the carrying capacity.
Taking these factors into account, the rate of cell proliferation in the continuum model is of the form, where S tot = n + r + e + m + p 1 b is the total cell density and p 1 b is the endothelial cell density corresponding to blood vessel surface density b, d cell is the carrying capacity for given cell type and H (d cell − S tot , ε d ) is a continuous approximation of the Heaviside function given as, is exactly equal to one, and equal to zero for (d cell − S tot ) ≤ −ε d . Also, the transition from 1 to 0 is smooth in the interval |d cell − S tot | < ε d .

Cell apoptosis
In the multiscale model, normal cell apoptosis depends on the cell's surroundings and intracellular levels of p53, denoted as [p53], updated at each simulation time-step, by integrating the following ordinary differential equations (ODEs) for time ∆t, with parameter values given in Table 1.
It must be noted that, following Alarcón et al (2005), k 8 is positive for normal cells and negative for cancer cells. This is because, in normal cells, p53 expression downregulates VEGF production. However, mutations of p53 in cancer cells leads to cases where p53 upregulates VEGF production (Alarcón et al, 2005, Royds et al, 1998 T HR , for a normal cell in presence of sufficiently many cancer cells, implying the assumption that the tumour microenvironment is hostile to normal cells (Alarcón et al, 2003, 2005, Gatenby and Gawlinski, 2001. Thus, the threshold [p53] T HR is given as, where ρ for a normal cell at site x is the proportion of normal cell density defined by, ρ = ∑ y∈θ x number of normal cells at site y ∑ y∈θ x number of normal cells + number of tumour cells at site y , where θ x is simply the cell's lattice site x if it contains more than one cell, and otherwise includes lattice sites from the Moore neighbourhood of x. In the continuum model we do not consider individual cells and their internal properties. Hence we model normal cell apoptosis as an oxygen-dependent rate, thus bypassing its dependence on the cell's internal p53 concentration. Solving Eq. (6a) at steady state, we get where α denotes high or low and subscript a indicates apoptosis. Substituting the parameter values given in Tables 1 gives, C high a = 1.5 mmHg, C low a = −7.4 mmHg.
As seen from Fig. 2a, Thus, in the absence of the tumour, normal cells undergo apoptosis when the oxygen concentration drops below C high a . Also, [p53] exceeds [p53] low T HR if C < C low a or C ≥ 0. Since C ≥ 0 in our model, in the tumour microenvironment, normal cells undergo apoptosis irrespective of the oxygen concentration.
Thus, assuming that the [p53]-[VEGF] subsystem in Eq. (6) is in quasi steady state, β norm is the rate of apoptosis, ρ n is the ratio of normal cell density to the sum of normal and tumour cell densities and following the above discussion, we conclude that the total rate of normal cell apoptosis in the continuum model is given by,   Note that since C high a denotes oxygen concentration for normal cell apoptosis, we will, henceforth, replace it by C norm a for notational convenience. Unlike normal cells, the tumour cells in the multiscale model are assumed to be more resistant to adverse conditions, like lack of oxygen. It was shown in Alarcón et al (2004) that the plot of division times of tumour cells against varying oxygen concentrations exhibits a vertical asymptote for a given oxygen concentration. This behaviour was not observed in the normal cells. The authors further state that such a behaviour of tumour cells is because of their ability of halting the cell cycle by entering the G0 phase rather that merely delaying or arresting it in the G1/S phase. The cell cycle model used in Owen et al (2011) is a simple phase model derived from the cell cycle model developed in Alarcón et al (2004). In Alarcón et al (2004), oxygen concentration affects the cell cycle via its regulation of p27 and cancer cells are assumed to enter or leave quiescence at certain threshold values of their p27 concentration. Since, Owen et al (2011) uses a simplified model, the p27 dependence is replaced with an equivalent dependence on the local oxygen concentration.
Thus, in the multiscale model, tumour cells enter into an intermediate state of quiescence (dormancy), when the oxygen concentration drops below C enter quiesc . During quiescence, they do not progress through the cell cycle and hence do not proliferate. If the oxygen supply improves and the concentration rises above C leave quiesc , they leave the state of quiescence and begin proliferating as usual. The quiescent tumour cells die if they remain in the state of quiescence for a duration longer than T quiesc .
In the continuum model, we account for quiescence by multiplying the proliferation rate with H (C − C q , ε q ), where C q is the mean of C enter quiesc and C leave quiesc . Additionally, we assume that the cancer cells undergo apoptosis at a rate given by H (C tum a −C)β tum where C tum a = C quiesc enter and β tum = 1/T quiesc . In summary, the equation for normal cell density is,

apoptosis in tumour microenvironment
where D norm is the random motility coefficient, T norm min is the minimum period of normal cell proliferation, C norm φ is the oxygen concentration at which the speed of normal cell proliferation is half maximal, d norm is the carrying capacity for cell proliferation and ρ n is as defined in Eq. (12).
The equation for tumour cell density is where D tum is the random motility coefficient, T tum min is the minimum time of tumour cell proliferation, C tum φ is the oxygen concentration at which the speed of tumour cell proliferation is half maximal and d tum is the carrying capacity for the tumour cell density.

Diffusibles
This layer accounts for diffusible molecules such as oxygen and the angiogenic growth factor VEGF. In the multiscale model, the diffusibles are modelled as continuous functions of space and time. As in Alarcón et al (2006) and subsequent articles, they are assumed to be in quasi steady state since the timescale for their diffusion on relevant tissue length-scale (minutes) is much shorter than the tumour doubling time (days/weeks). The general form of the equation describing the evolution of a diffusible, U, given in Owen et al (2011), is where D u is the diffusion coefficient of U in the extracellular space, ψ u (x,t) is the vessel permeability to U, U blood is the concentration of U in the blood, S u (x,t) denotes the cell and environment dependent production/removal rate of U, δ u is the rate of decay of U and ρ u (x,t) denotes the vascular surface density at x defined as the surface area of the cylindrical vessel of radius R(x,t) and length L(x,t), divided by the volume of the lattice site with sides ∆ x. If no vessel is present at a lattice site, ρ u is set to zero. A finite difference approximation, on the same lattice as the cells, is used to approximate solutions to Eq. (15). For oxygen, in the multiscale model, vessels act as sources and deliver oxygen to the surrounding tissue. Cells consume oxygen thus acting as sinks. Also, oxygen does not decay naturally implying δ c = 0. If two or more cells are present at the same lattice site, each cell makes a contribution to S c (x,t). The uptake/secretion rates are given as rates per cell.
For VEGF in the multiscale model, the level of VEGF in the blood is assumed to be negligible, V blood = 0. Also, VEGF decays naturally at a rate given by δ v . Normal cells secrete VEGF when the intracellular VEGF level, [VEGF], updated by solving the ODEs in Eq. (6), exceeds the threshold [V EGF] T HR . Since intracellular processes are not included in the continuum model, we replace the [VEGF] dependence with an equivalent dependence on local oxygen concentration. As seen from Eq. (6b), oxygen promotes degradation of intracellular VEGF. Hence, as done for normal cell apoptosis, solving the [p53]-[VEGF] subsystem in Eq. (6) at steady state, for C, with [V EGF] = [V EGF] T HR and other parameter values as given in Table 1, we get C = 1.7 and 3.8. It can be also seen from Fig. 2, that for C between 1.7 and 3.8, [V EGF] > [V EGF] T HR at steady state. Since cells die when the oxygen concentration drops below the lower threshold, 1.7 mmHg, we consider C norm v to be equal to only the upper threshold of 3.8 mmHg.
Additionally, it is assumed in the multiscale model that tumour cells secrete VEGF when in quiescent state. Since, in the continuum model, C q denotes the oxygen concentration below which tumour cells enter quiescence, we use the same threshold to model VEGF secretion by tumour cells.
Thus, in the continuum model, we model the diffusibles as Eq. (15). The uptake and secretion rates per cell from the multiscale model are converted to rates per unit cell density.
The equation for oxygen concentration is therefore, where D c is the diffusion coefficient, P c is the permeability of vessels to oxygen, C blood is the oxygen tension in the blood and k cell c is the rate of oxygen consumption by cell of type cell.
The equation for VEGF concentration is, where D v is the diffusion coefficient, P v is the permeability of vessels to VEGF, k cell v is the rate of VEGF secretion by cell of type cell and δ v is the rate of VEGF decay.

Vasculature
The vascular layer represents the blood vessel network. In the multiscale model, at each time step, each vessel segment undergoes structural adaptation: changes in the vessel radius, blood flow and pressure, and vessel pruning in response to stimuli including haemodynamic and metabolic effects. A detailed description of the stimuli and the process of adaptation can be found in Owen et al (2009Owen et al ( , 2011. It is far from clear how to include such local dynamics in a continuum model. Hence, for the first approximation, we assume fixed vessel radii in all the simulations for simplicity and ease of comparison. Additionally, this layer also accounts for angiogenesis, the process of formation of new blood vessels. In the multiscale model, endothelial tip cells sprout from preexisting blood vessels, in a time step of length ∆t, with a probability where x denotes the lattice site at which the vessel cell is currently located, P max sprout is the maximum probability of sprouting per unit time, and V sprout is the VEGF concentration at which the probability is half maximal. Also, tip cells cannot sprout if the number of cells at x exceed the carrying capacity d tip .
Thus in the continuum model, we model the rate of tip-cell generation as Tip cell movement, in the continuum model is given as chemotactic response in the direction of increasing VEGF which is similar to Eq. (2), the continuum limit of Eq. (1), the probability of cell movement in the multiscale model. We multiply the chemotaxis flux term by S , which is zero when S tot > d tip and one otherwise, to ensure that, in the numerical simulations, the flux is zero when S tot > d tip .
Since in the multiscale model, the tip cells can form connections only after having moved a certain distance away from the parent blood vessel, we assume that in the continuum model, the tip cells, e(x,t), sprouting from the pre-existing vessels are immature and convert into mature tip cells, m(x,t), at a constant rate. These mature tip cells are capable of forming new connections with other mature tip cells and preexisting vessels. The tip kinetics of forming new vessels are modelled as where the first term describes the formation of loops created by two capillary tips (tip-tip anastomosis) at a rate λ 1 , while the second term describes the joining of a tip to the side of a capillary (tip-capillary anastomosis) at a rate λ 2 (Aubert et al, 2011, Balding and McElwain, 1985, Edelstein-Keshet and Ermentrout, 1989, Gaffney et al, 2002.
Without structural adaptation and vessel pruning, and assuming fixed vessel radii, the blood vessel density, in the multiscale model, changes only as a result of new connections formed through angiogenesis.
Thus, in the continuum model, the evolution equations for immature tip cell density e, mature tip cell density m, and vessel density b are given by tip-tip and tip-capillary anastomosis , where D e is the diffusion coefficient, χ e is the chemotactic coefficient, d tip is the carrying capacity for tip cell movement and sprouting, α 1 is the rate of conversion of immature tip cells to mature ones, p 1 is endothelial cells per unit vascular surface area used as a conversion factor to calculate endothelial cell density corresponding to a given blood vessel surface density. A detailed discussion of the parameter values can be found in Appendix A. Some of the parameter values, for instance the diffusion coefficients and permeability of VEGF, were taken directly from Owen et al (2011). Others, like carrying capacity with units of number of cells, and oxygen consumption rate and VEGF secretion rate with units of per cell per unit time, from Owen et al (2011), had to be converted into suitable units, like number of cells per unit volume or per unit cell density per unit time, for use in the continuum model.
The multiscale model is implemented in C++ using CVODE to integrate subcellular ODEs, and SuperLU to solve linear systems for flow calculations and reactiondiffusion equations. Other computational details can be found in the Supplementary material of Owen et al (2011). To solve the system of PDEs representing the continuum model, we use the D03PCF routine from the NAG library which employs the method of lines to reduce the PDEs to a system of ODEs and solve the resulting system using a backward differentiation formula method.

Simulation of a normal vascular tissue
In the absence of tumour cells, the continuum model equations reduce to Eqs. (32) given in Appendix B.
The multiscale model has been implemented in 2D  and 3D (Perfahl et al, 2011). However, in this paper, we establish scenarios that enable comparison with the 1D continuum model. Therefore, we consider a rectangular domain [0, L x ] × [0, L y ] with L x > L y and take averages of the simulation results in the ydirection. Fig. 3 shows the initial tissue configurations from the two models. The multiscale model consists of one blood vessel situated close to the left boundary and surrounded by only normal cells. No tip cells are present. Oxygen and VEGF are initialised to their respective quasi-steady states. The cell and vessel densities computed from the multiscale model by averaging in the y-direction are used to set corresponding initial conditions for the continuum model. Oxygen and VEGF concentrations are then calculated from Eqs. (32a), (32b) respectively. In both models, the oxygen concentration is higher in regions of high vessel density owing to delivery of oxygen from the vessels. Normal cells secrete VEGF in the hypoxic regions. VEGF promotes sprouting of the tip cells. The tip cells form new vessels and the network expands into the region with low vessel density.

Vessels
As seen from Fig. 5, the density profiles for all the variables are comparable in the two models. However, at t=3 days, in the continuum model, the vessel density rises only at the site of anastomosis (the small peak seen on the left of the pre-existing vessel). However, the multiscale model uses a snail-trail approach wherein the endothelial cells left behind by the moving tip cells also contribute to the vessel density on anastomosis. Therefore, we see a continual change in the vessel density, starting from the pre-existing vessel towards x=0. Also, at these sites of anastomosis, the oxygen concentration rises above C norm v , the threshold for VEGF secretion by normal cells. Hence, tip cell sprouting in this region reduces, as seen from Fig. 5. The tip cells already present in the tissue move chemotactically in the direction of hypoxia (up the VEGF gradient). Therefore, at t=5 days, in the continuum model the vessel density to the left of the pre-existing vessel remains almost unchanged. However, in the multiscale model, we notice a small rise owing to the endothelial cells forming new vessels. The continuum model seems to agree well with the multiscale model, except for the discrepancies seen in vessel density and consequently in oxygen concentration. We believe that the snail-trail approach from the multiscale model is not accounted for in the continuum model and might be the cause of these differences. 4 Steady states of the spatially homogeneous normal tissue submodel As seen from Fig. 4, the normal cell density, in the continuum model, exhibits a travelling wave moving from the healthy and well vascularised part of the tissue to a part of the tissue without any cells and blood vessels. Since travelling wave solutions connect two spatially homogeneous steady states, we determined such steady states of the system given by Eqs. (32). The analysis also enables us to understand which of these steady states connect to form the waves seen in the simulations and investigate the wave speed and its dependence on the model parameters and initial conditions of model variables.
Based on the analysis detailed in Appendix C, we list below four types of spatially homogeneous steady state of Eqs. (32).
This trivial steady state represents a tissue with no cells, no blood vessels and an arbitrary oxygen concentration C * ≥ 0.
Another trivial steady state representing a tissue with only blood vessels of arbitrary density, b * > 0 and oxygen supplied and consumed by those vessels.
A non-trivial steady state representing a tissue with normal cells, blood vessels and oxygen. The sum of normal cell density and vessel density exceeds the carrying capacity for normal cells thus blocking normal cell proliferation. Also the vessel density is such that it supplies enough oxygen (C > C norm v + ε v > C norm a + ε a ), so that there is no normal cell death, no VEGF production, and hence no angiogenesis. 4. m = 0, e = 0, n = n * , b = b * such that n * + (p 1 b * ) > d tip + ε d , This steady state, also non-trivial, represents a tissue with normal cells, blood vessels, oxygen and VEGF but no tip cells. The oxygen attains a certain concentration dictated by the blood vessel density and the normal cell density. This value of the concentration being below the threshold C norm v + ε v , causes hypoxic normal cells to secrete VEGF. However, the normal cells and the blood vessels occupy the entire space in the tissue thus preventing any further sprouting from the blood vessels in response to the VEGF signals. Similar to the oxygen concentration, the VEGF concentration reaches a certain steady concentration depending on the normal cell density and the blood vessel density.
The non-trivial steady states suggest that n and b span a family of steady states such that n + (p 1 b) > d norm + ε d . The steady state oxygen concentration, given by Eq. (44b), is a function of the normal cell density and the blood vessel density. Tip cell densities are zero at steady state and VEGF concentration is also zero in all cases, except when n + (p 1 b) > d tip + ε d and C norm a + ε a < C < C norm v + ε v . To illustrate how different initial conditions lead to different steady states, we simulated the system of DAEs given by Eqs. (33) in Appendix C, using the MATLAB routine ode15s, with different initial values of n and b and for the parameter values given in Appendix A. The time-dependent solutions for each of these initial configurations projected on the (n, b) plane are shown in Fig. 6. As seen from Fig. 6, the (n, b) space can be divided in eight different regions depending on the thresholds of oxygen for normal cell apoptosis and VEGF secretion and of total cell density for normal cell proliferation and tip cell production. In Fig.  7, we show the temporal evolution of all variables for trajectory a, which traverses the largest number of regions in Fig. 6.
We will now discuss how the path of the trajectory and thence the steady state configuration attained by the system varies for initial conditions lying in each of these regions.

I Apoptosis, VEGF secretion -S tot
The initial VEGF concentration is non-zero, since C < C norm a + ε a < C norm v + ε v . However, since the total cell density exceeds the carrying capacity of tip cell sprouting, no new tips are formed. Normal cells undergo apoptosis due to hypoxia. Trajectories from this region (e.g: trajectories a and b) enter region II. II Apoptosis, VEGF secretion, angiogenesisd norm + ε d < S tot < d tip + ε d , C < C norm a + ε a In this region too, since C < C norm a + ε a initially, the normal cell density drops due to cell apoptosis. Tip cells are formed in response to VEGF. However, normal cells do not proliferate since the total cell density exceeds its carrying capacity. Depending on the vessel density, the trajectory enters region III (trajectory a) or V (trajectory b). III Apoptosis, proliferation, VEGF secretion, angiogenesis -S tot < d norm + ε d , C < C norm a + ε a Following trajectory a through this region, we see that the normal cell density continues to drop due to hypoxia. However, since the total cell density is less than the carrying capacity for tip cell sprouting, new tip cells form in response to VEGF. This leads to an increase in the blood vessel density followed by an increase in the oxygen concentration. This rise in oxygen and availability of space, since S tot < d norm + ε d , promotes normal cell proliferation. Once the oxygen concentration rises beyond C norm a + ε a , apoptosis stops and the trajectory enters region IV. IV Proliferation, VEGF secretion, angiogenesis - vessel density and thus oxygen concentration within the tissue. This rise in oxygen promotes further normal cell proliferation until S tot crosses d norm + ε d . Trajectories from this region (e.g.: trajectories a and d) enter region V V VEGF secretion, angigoenesisd norm + ε d < S tot < d tip + ε d , C norm a + ε a < C < C norm v + ε v As seen from trajectory c and the continuation of trajectory a in this region, normal cells do not proliferate as the total density exceeds their carrying ca- pacity. Since C < C norm v + ε v , VEGF is secreted by the normal cells and new tip cells are formed leading to a further rise in blood vessel density and the oxygen concentration. This process continues until the oxygen concentration exceeds C norm v + ε v , VEGF secretion stops, and the trajectory enters region VIII. VI VEGF secretion -S tot > d tip + ε d , C norm a + ε a < C < C norm v + ε v For initial conditions in this region (e.g.: trajectory g), there is enough oxygen for normal cells to survive since C > C norm a + ε a . However, they do not proliferate due to lack of space. Thus the normal cell density remains unchanged. Normal cells secrete VEGF since C < C norm v + ε v , but no new tip cells can be formed due to unavailability of space. Hence new blood vessels cannot be formed and the blood vessel density also remains unchanged. All variables are in a steady state with non-zero VEGF concentration.

VII Proliferation -S tot
+ ε a , normal cells do not secrete VEGF and do not undergo apoptosis. No VEGF implies the vessel density remains unchanged. However, normal cells proliferate and hence their cell density increases until S tot crosses d norm + ε d , proliferation stops, and the trajectory (e.g.: trajectory e) enters region VIII. However, in cases like trajectory d with low blood vessel density and hence low oxygen concentration, the increase in consumption of oxygen, due to the increase in normal cell density, causes the oxygen concentration to drop below the threshold C norm v + ε v . In this case, the trajectory traverses through regions IV and V before entering region VIII. VIII Steady state -d norm + ε d < S tot < d tip + ε d , C > C norm v + ε v In this region, normal cells do not proliferate due to lack of space. No VEGF is secreted by the normal cells since C > C norm v + ε v . Hence, angiogenic sprouting does not occur. For initial conditions in this region (e.g.: trajectory f) the variables are already in the steady state and hence the trajectory does not move at all. For trajectories entering from region V or VII (e.g: trajectories a to e), the tip cells already present within the tissue continue to mature and form connections via anastomosis. Hence the blood vessel density continues to rise until all the tip cells form connections. Once this process ends, all the variables attain a steady state and the trajectory stops. IX Steady state -S tot > d tip + ε d , C > C norm v + ε v For initial conditions starting in this region (e.g.: trajectory h), the variables are already in steady state since cells do not proliferate due to lack of space and also do not secrete VEGF or promote angiogenesis due to sufficient availability of oxygen. No trajectory starting in any other region can enter this region since it will eventually enter region VIII and attain a steady state as described earlier.
Since any initial configuration of n and b would lie in one of the nine regions mentioned above, the course of the solution will qualitatively follow the course of one of the trajectories shown in Fig. 6. However, it must be noted that the simulations in this section were performed with parameter values corresponding to the published multiscale model simulations as given in Appendix A. Changing the parameter values changes the behaviour of the system. For instance, as seen from trajectory a, we see an instantaneous drop in the normal cell density. This is because the rate of normal cell death (β norm = 1 min −1 ) is much higher than the rate of tip cell sprouting P max sprout V V sprout +V = 1.67e − 4 min −1 , with V = 1 Nm, the typical VEGF concentration assumed in the model. For β norm < P max sprout V V sprout +V trajectory a follows an altered path shown by the red dashed line in Fig. 6. In this case, because tip cells and eventually blood vessels are formed at a rate faster than normal cell apoptosis, the oxygen concentration rises beyond C norm a + ε a before trajectory a enters region III. As a result, trajectory a continues from region II into region V.

Wavespeed of normal tissue expansion for different initial conditions
We now seek for an estimate of the wavespeed of a travelling wave solution to Eqs. (32) given in Appendix B. As x → ∞, we have n → 0 b → b ss , C → C ss , V → 0, e → 0 and m → 0. Linearising Eq. (32c) around this state, we obtain Eq. (22) is qualitatively similar to the Fisher-KPP equation, linearised at the leading edge, given as, where α is the growth rate, K is the carrying capacity and D is the diffusion coefficient.
Previous studies (Fisher, 1937, Larson, 1978, McKean, 1975, Mollison, 1977 have shown that Eq. (23) admits a travelling wave solution of the form u(x,t) = f (x − ct), for some function f , where c denotes the speed of the wave. Additionally, for initial conditions satisfying u → K as x → −∞ and u = O(e −µx ) as x → ∞ the travelling wave evolves with a speed given as By analogy with Eq. (24), we expect that a travelling wave whose normal cell density profile decays as e −µx as x → ∞ has the approximate wavespeed To test how the wavespeed for normal cell density varies with µ, the rate at which n → 0, we ran simulations with initial condition of normal cell density given as where x 0 = 1 mm. The blood vessel density was set to 20 mm 2 /mm 3 homogeneously, tip (immature and mature) cell density was zero. Solving Eq. (32a) and Eq. (32b), at the leading edge of the wave where n ss → 0, e ss = m ss = 0, and b ss = 20, we get Also note that at the leading edge, H (d norm − S ss tot ) = 1 since S ss tot = p 1 b ss < d norm . Thus, by substituting C from Eq. (27a) in Eq. (25c), we can determine a µ , a min and µ 0 , the value of µ such that the wavespeed is minimum for µ ≥ µ 0 . Black circle: a µ given by Eq. (25a), the speed associated with the initial condition decay rate µ. Blue cross: a min , the predicted minimum speed given by Eq. (25b). Red cross: simulated wave speed.
As seen from the Fig. 8, for smaller values of µ, the simulated speed equals a µ , given by Eq. (25a). However, as µ exceeds α norm (C ss ) D norm = 38.7 for C ss obtained from Eq. (27a), the normal cell density wave attains the minimum speed given by Eq. (25b). Note that C ss given by Eq. (27a) remains constant when e ss = m ss = 0, n ss → 0 and b ss is any arbitrary non-zero value less than d norm . Consequently, a norm is constant for such a setup implying that the speed for normal cell density wave depends only on µ.
If the initial vessel density is such that p 1 b ss ≥ d norm , then C ss is still given by Eq. (27a). However, a norm will be zero since H (d norm − S ss tot , ε d ) = 0 and H (C norm a − C, ε a ) = 0. Thus, the normal cell density wave does not move in such a case. In our (a) Top: Initial normal tissue in the multiscale model consisting of one blood vessel surrounded by normal cells. Bottom: Tumour cells introduced in the simulated well vascularised normal tissue. This configuration is used as initial condition for tumour growth simulation. Note that we ran 10 such realisations of a well vascularised normal tissue, which were then used as initial configurations for 10 realisations of the tumour growth simulation.  model simulations, we observe only a small change in the initial normal cell density profile owing to diffusion by random motility (result not shown).

Tumour growth within a normal vascular tissue
In this section, we perform simulations of a tumour tissue. Henceforth we will use the term continuum model for the full system given by Eqs. To study tumour growth, we first simulate a well vascularised normal tissue, in the multiscale model, starting with the same setup described in Sect.
(3). The initial and the final configuration (850 days) of the lattice is presented in Fig. 9a. Tumour cells are introduced to the final configuration and the simulation is continued, similar to simulations done in Owen et al (2009Owen et al ( , 2011. In the continuum model, we set the normal cell density, homogeneously, equal to its carrying capacity, d norm , and the vessel density, homogeneously, to the average vessel density calculated from the multiscale model at the end of 850 days. No tip cells, mature or immature are present at the start of the simulation. All the parameter values are as defined in Appendix A. The density of tumour cells obtained from the multiscale model, by averaging in the y-direction, is used as initial condition in the continuum model. The initial densities of tumour cells, normal cells and the blood vessels, in both the models, are presented in Fig. 9b.

Multiscale Model
Continuum Model  10 shows the spatio-temporal evolution of the normal and tumour cells. Tumour cells create a microenvironment that is hostile to normal cells and triggers their apoptosis. This creates space for tumour cell proliferation and invasion of the healthy tissue. In this particular simulation setup, the tissue is already sufficiently well vascularised to supply enough nutrients for tumour growth. Hence no new vessels are formed.

Wavespeed of tumour invasion for different initial conditions
The linearised form of Eq. (14) about the leading edge of the tumour invasion wave, where r → 0, n → n ss , b → b ss , C → C ss , V → 0, e → 0, m → 0, (n ss , b ss and C ss are values attained by n, b and C as x → ∞), is given by Thus, as done for normal tissue expansion, the wave speed for tumour cell density with initial conditions that decay as e −µx as x → ∞ is given by To test how the initial condition decay rate, µ, affects the wave speed, we ran simulations for tumour invasion of a normal vascular tissue with initial tumour cell density given as where x 0 = 1 mm. Additionally, n(x, 0) = d norm , b(x, 0) = 20 mm 2 /mm 3 , e(x, 0) = m(x, 0) = 0. Similar to normal tissue expansion, the oxygen and VEGF concentration at the leading edge of the tumour wave are given by C ss = P c C blood b ss b ss (k vessel c p 1 + P c ) + k norm c n ss (31a) Also, H (d tum − S ss tot , ε d ) = 1 since S ss tot < d tum . We can hence determine a tum , a tum µ , a tum min and µ 0 , the value of µ after which the tumour wave attains the minimum speed. Fig. 11 shows the tumour cell density wave speed for different values of µ. The simulated speed of the wave decreases with increasing µ. For µ smaller than the threshold α tum (C) D tum , determined to be 35.2 for C ss given by Eq. (31a), the wave speed equals a tum µ given by Eq. (29a) and attains the predicted minimum speed, a tum min given by Eq. (29b) as µ exceeds the threshold. Wave Speed (mm/day) Fig. 11: Wave speed for the tumour cell density, obtained from the simulation of the continuum model for initial condition given in Eq. (30) for different values of initial condition decay rate µ. The simulated speed decreases with increasing µ. The simulated speed equals a tum µ until it crosses the threshold α tum (C) D tum , determined for C ss given by Eq. (27a). For µ higher than the threshold, the simulated speed attains the predicted minimum speed, a tum min . Black circles: a tum µ , given by Eq. (29a), the speed associated with initial conditions decay rate, µ. Blue crosses: a tum min given by Eq. (29b) is the predicted minimum speed. Red crosses: simulated tumour cell density wave speed.
Unlike the case of normal tissue expansion, C ss , given by Eq. (31a), depends on the normal cell and vessel densities at the leading edge of the tumour cell density wave. Thus, the wave speed for tumour cell density depends not just on µ, the initial condition decay rate for tumour cell density, but also on the normal cell and vessel densities at the leading edge of the tumour wave. Thus, the speed of invasion of the normal tissue by the tumour cells depends on the configuration of the normal tissue.
6.2 Dependence of wavespeed on the carrying capacity of tumour cell density As mentioned in Sect. (6.1), the oxygen concentration at the leading edge of the tumour cell density wave does not depend on the densities of other model variables. Hence, the wave speed for tumour cell density given by Eqs. (29), depends only on the decay rate, µ, and the availability of space given by H (d tum − S ss tot , ε d ). For S ss tot ≥ d tum +ε d the minimum wave speed predicted by Eq. (29b) is zero. To determine the dependence of the wavespeed on the carrying capacity, we ran simulations of both the continuum and the multiscale models for different values of d tum . The initial conditions in the continuum model were as described in Sect. (6.1). We choose µ = 10 5 since it approximates to the step function initial conditions used in the multiscale model and the wave associated with it attains the minimum speed.
As seen from Fig. 12, the wave speed decreases as the carrying capacity of the tumour cell density decreases (i.e. as d tum − S ss tot decreases). Note that, for d tum ≤ S ss tot + ε d , the simulated speed is not equal to the predicted minimum speed of zero. . Green dashed line: S ss tot = d tum + ε d , the point of transition from pushed to pulled waves. For d tum < S ss tot − ε d , the tumour wave is a pushed wave since the speed of the tumour cell density wave depends on the whole wave front rather than the leading edge of the wave. Consequently, the (minimum) simulated speed of the wave is greater than the predicted minimum speed of zero. Stokes (1976) refers to such waves as pushed waves since their speed of propagation is determined by the whole wave front and not by the leading edge, as in the case of pulled waves. In our model simulation, the tumour cells at the wave front create a hostile environment for the normal cells there, thus triggering their apoptosis. As the normal cells die, S tot drops below d tum + ε d consequently enabling tumour cell proliferation and progression of the wave. Thus, when S tot > d tum + ε d , the speed of the wave for the tumour cell density does not depend on the leading edge of the wave but on the whole wave front of the tumour cell density.

Discussion
The focus of this work has been to develop a continuum model based on the multiscale model given in Owen et al (2011) for growth of a vascular tumour tissue. Both the models account for evolution of normal and tumour cells and formation of vessels in response to surrounding conditions like availability of oxygen and secretion of VEGF. The results from the 1D continuum model were compared to the averaged density profiles obtained from a long 2D strip in the multiscale model. As seen from the density profiles in Fig. 4 and Fig. 10 the model results are qualitatively similar. However, in Fig. 4, the vessel density away from the parent vessel is low ( 5 mm 2 /mm 3 ) in the continuum model compared to that in the multiscale model ( 20 mm 2 /mm 3 ). One plausible reason for the difference in the vessel density could be the use of the snail-trail approach for movement of tip cells in the multiscale model. In this approach, it is assumed that the tip cells leave a trail of endothelial cells while moving within the tissue. On anastomosis, the entire trail forms a new vessel thus adding to the vessel density. We believe that we do not account for the snail trail approach in continuum model. However, we can extend the model to account for the snail trail approach as done in Balding and McElwain (1985), Byrne and Chaplain (1995), Connor et al (2015), Spill et al (2015).
Moreover, not all processes from the multiscale models can be translated into rates in the continuum model. For instance, in the multiscale model, normal cell death and VEGF secretion depend on the intracellular p53 and VEGF concentration which evolve over time according to the local oxygen supply. We bypassed this subcellular dependence in the continuum model by assuming that the p53-VEGF subsystem attains a quasi steady state, so the thresholds in terms of p53 and VEGF become thresholds in terms of local oxygen concentration. A parameter sensitivity analysis, where parameters in the continuum model, not directly taken from the multiscale model, are varied within physiological ranges, could provide additional information about the relative impact of individual parameters on the continuum model dynamics.
A drawback of the continuum model is that it is not clear how to include certain local dynamics like blood flow and adaptation of vessel radius to the microenvironment. However, with the continuum model, it was possible to analyse the steady states of the system, learn about wave speeds and their dependence on the model parameters and initial conditions. We determined that n and b span a family of steady states and showed that the steady state attained by the system depends on initial conditions of the model variables. As seen from Fig. 6, the n − b plane is divided into nine distinct regions by the thresholds of oxygen for normal cell apoptosis and VEGF secretion and of total cell density for normal cell proliferation and tip cell sprouting. Any initial condition for n and b will lie in one of these regions and thus qualitatively follow the course of the solutions shown in Fig. 6.
The system given by Eqs. (13)- (14), (16)- (17) and (19)-(21) admits travelling waves for normal tissue expansion and tumour invasion into normal tissue. In the case of normal tissue expansion, the oxygen concentration at the leading edge of the wave remains constant for any non zero value of the blood vessel density. Thus, the speed of the wave depends only on the decay rate of the initial profile of normal cell density. The wave for tumour invasion also depends on decay rate of the initial profile of tumour cell density. Additionally, it can be seen from Eq. (31a) that the oxygen concentration at the leading edge of the wave depends on the normal cell and vessel densities at there. Thus, the speed at which the tumour invades a normal vascularised tissue depends on configuration of that tissue. We also showed that the speed of the invading tumour wave varies with d tum , the threshold of total cell density for tumour cell proliferation and changes from a pulled to a pushed wave if d tum is less than the total cell density at the leading edge of the wave. These pushed waves of tumour invasion propagate by creating a micro-environment that is hostile for the normal cell survival, thus triggering apoptosis of normal cells at the wave front.
To the best of our knowledge, a model as complex as the multiscale model from Owen et al (2011) has not been compared to an equivalent continuum model before. We remark that the focus of this paper was on the published multiscale model from Owen et al (2011). However, it would be interesting to consider different models for cell cycles, alternative rules for cell apoptosis and entry to/exit from quiescence, in the multiscale model and then correspondingly in the continuum model.
The multiscale model from Owen et al (2011) also includes conventional drug therapy and a macrophage based gene therapy. In the multiscale model, the chemotherapeutic drug is assumed to be a continuous variable and its evolution is given by Eq. (15). It would thus be straightforward to include a chemotherapeutic drug as an additional variable in the continuum model. Cells, in the multiscale model, with local drug concentration higher than a prescribed threshold intercalate the drug and die on attempting to divide. This effect of the drug on the cell densities in the continuum model can be modelled using H given by Eq. (5). For the macrophage based gene therapy in the multiscale model, macrophages are assumed to extravasate from the pre-existing blood vessels at a certain probability, similar to tip cell extravasation. Hence, as done for tip cell sprouting, we can readily model macrophage extravasation in the continuum model. The extended continuum model can then be used to analyse the impact of drug and macrophages on the speed of the tumour invasion wave. Such analysis might further help in designing novel therapeutic strategies to reduce the speed at which the tumour invades a normal vascularised tissue.

A Parameter values
In the multiscale model by Owen et al (2011), a lattice site corresponds to a cube with length of each side equal to ∆ x = 40µm = 0.04 mm. Let us denote the volume of one lattice site by V site = 6.4 × 10 −5 mm 3 . For the remainder of this section, lattice site has the same meaning as above unless stated otherwise.

A.1 Estimation of P c
We ran simulations of the multiscale model for different vessel densities on a 1 × 1 square lattice with only one normal cell. The simulations were run only for one time step and the oxygen concentration was noted. Since the lattice dimensions are 1 × 1, we can assume that the system is spatially homogeneous and thus obtain P c by substituting the values of n, b and C in Eq. (16). The value is given in Table4 Table 3 A

.3 Estimation of d cell
In the multiscale model, the carrying capacity of a cell at lattice site x is defined as the maximum number of cells that can be present at x. Therefore, we obtain the carrying capacity for the continuum model as number of cells per unit volume (cells/mm 3 ) by dividing the carrying capacity by V site . The value is given in Table 3.

A.4 Estimation of D cell
In the multiscale model, apart from random movement, normal and tumour cells also spread by via cell division by placing their daughter cells in neighbouring sites if the number of cells in the parent cell's site exceeds its carrying capacity. A cell of type cell takes at least T cell min minutes to divide. On division, the new cell is placed in one of Moore neighbourhood sites. Thus, the mean squared displacement of the cell in T cell min min is ∆ x mm. For normal cells, the random motility coefficient, in the multiscale model is zero. We thus estimate the value of D norm using the mean squared displacement estimated above. For tumour cells, the random motility coefficient, in the multiscale model is non-zero. Thus, starting with the estimate obtained from the mean squared displacement of tumour cells, we simulate the continuum model for different values of D tum and compare the distance covered in unit time with the simulation from the multiscale model. The best fit thus obtained was chosen as the value for D tum .
All the values of D cell are given in Table 2. A.5 Estimation of p 1 In all simulations of the multiscale model presented in this paper, all vessels have a fixed radius R = 0.006 mm and length L = ∆ x. Also, in the multiscale model, all blood vessels in a lattice site are associated with one endothelial cell. Therefore in a lattice site with volume V site , each vessel with surface density of 2πR∆ x/V site = 23.56 mm 2 /mm 3 has an endothelial cell density of 1/V site = 15625 cells/mm 3 associated with it. Therefore, p 1 defined as the number of endothelial cells per unit vascular surface area is obtained by dividing the endothelial cell density by the vessel surface density. Thus p 1 = 663 cells per mm 2 , as given in Table 5. A.6 Estimation of λ 2 Estimating the rates of tip-tip and tip-vessel anastomosis from the multiscale model is difficult. We thus use the values given in Schnepf et al (2008). To convert λ 2 from units of per unit capillary length density per unit time to the units of per unit capillary surface density per unit time, we multiply λ 2 from Schnepf et al (2008) with the surface area of a vessel with radius 0.006 mm and unit length. The values of λ 1 and λ 2 are given in Table 5.  (14), (16)- (17) and (19)- (21), that can be taken directly from the multiscale model  and are related to cell proliferation (sprouting for tip cells), apoptosis, movement, oxygen consumption and VEGF secretion. Dash indicates that the parameter is not defined for that cell type.   (14), (16)- (17) and (19)- (21), that are modified from those of the multiscale model  and are related to cell proliferation (sprouting for tip    Rate of tip-tip anastomosis (Schnepf et al, 2008) λ 2 (mm 2 /mm 3 ) −1 min −1 2.8 × Blood Vessel Density ∂ b ∂t = 1 p 1 (λ 1 m 2 + λ 2 mb) tip-tip and tip-capillary anastomosis . (32f)

C Steady state analysis
Assuming all the variables to be spatially homogeneous, we obtain the following set of differential algebraic equations (DAEs), 0 = P c b(C blood −C) implying e = 0. Substituting e = 0 in Eq. (33d) gives Note: In the remainder of this section, S tot = n + (p 1 b), since e = 0 and m = 0. Eq. (35) holds if, either where H (z, ε) is as defined in Eq. (5). The smoothness of transition from 0 to 1 is given by the value of ε. As the value of the function is zero for z ≤ −ε, steady state will only be attained when z ≤ −ε. For ε = 0, the function reduces to the discontinuous Heaviside function reflecting the dynamics observed in the multiscale model. We now analyse the possible cases leading to steady states of Eqs.
This steady state represents a tissue with no cells and no blood vessels and an arbitrary oxygen concentration.
This state is thus a special case of Case 1.1 representing a tissue with no cells, no blood vessels and no oxygen.
Case 2 When V = 0, Eq. (33b) reduces to It must be noted that we have assumed V blood = 0. This leads to two possibilities, either C > C norm v + ε v or n = 0. Case 2.1 If n = 0, Eq. (33a) becomes which implies C = P c C blood k vessel c p 1 + P c . Also n = 0 satisfies Eq. (33c). Thus the next possible steady state configuration is given by m = 0, e = 0, V = 0, n = 0, b = b * > 0, C = P c C blood k vessel c p 1 + P c .
This represents a tissue with only blood vessels of arbitrary density, b * > 0 and oxygen supplied and consumed by those vessels.

Case 2.2
Let us now consider the other possibility, C > C norm v + ε v . Since C norm v + ε v > C norm a + ε a , we have C > C norm v + ε v > C norm a + ε a . Hence, Eq. (33c) reduces to ln(2) T norm For Eq. (42) to be true, either n = 0 or S tot > d norm + ε d . Case 2.2.1 Let n = 0. This leads us back to Case 2.1, implying that the steady state configuration is given by Eq. (41).
Although for the given parameter set, C = P c C blood k vessel c p 1 + P c > C norm v + ε v holds, the inequality is irrelevant in this case because, irrespective of the oxygen concentration, there will be no VEGF production in the absence of normal cells (n = 0). Case 2.2.2 Next, let us assume, S tot = n + (p 1 b) > d norm + ε d . From Eq. (33a) we get, n + (p 1 b) > d norm + ε d , where I(X) = P c C blood − X k vessel c p 1 − X P c X k norm c . ( This state represents a tissue with normal cells, blood vessels and oxygen. The sum of normal cell density and vessel density exceeds the carrying capacity for normal cells thus blocking normal cell proliferation. Also the vessel density is such that it supplies enough oxygen (C > C norm v + ε v > C norm a + ε a ), so that there is no normal cell death and no VEGF production leading to angiogenesis.