ࡱ> 68J612345_dfhjl n p r t y{}ACEGIKMOQSUWY[] _!!a""c##e$$g%%i&&k''m((o))q**s++u,,w--y..{//}0012A22C33E44G55#` 'bjbj P ?``` j 222F8jLFif j,,,T4DFFFFFF$ϴh7\j2<+^<+<+jvv,TK#@d@d@d<+2#v,(2TD@d<+D@d@dRT62\T  00'(k nNF H 90ifZ" <\2\D }#@d% 'TQ,jjcj}}}i<+<+<+<+FFF`FFFFFFvvvvvv Multiscale Modelling and Nonlinear Simulation of Vascular Tumour Growth Paul Macklina, Stephen McDougall*b, Alexander R.A. Andersonc, Mark A.J. Chaplainc, Vittorio Cristinid, John Lowengrub*a,e aDepartment of Mathematics, University of California Irvine, Irvine, CA 92697-3875, USA bDepartment of Petroleum Engineering, Heriot-Watt University, Edinburgh EH14 4AS, Scotland cDivision of Mathematics, University of Dundee, Dundee DD1 4HN, Scotland dBiomedical Engineering, School of Health Information Sciences, University of Texas Health Science Center, Houston, TX 77225, USA eDepartment of Biomedical Engineering, University of California Irvine, Irvine, CA 92697 Abstract The growth and development of solid tumours occurs in two distinct stages - the avascular growth phase and the vascular growth phase. During the former, the tumour remains in a diffusion-limited, dormant state of a few millimetres in diameter (cf. multicell spheroids, carcinoma in situ) while during the latter, invasion and metastasis may take place. In order to accomplish the transition from avascular to vascular growth, solid tumours secrete diffusible substances known as tumour angiogenic factors into the surrounding tissue. These factors in turn induce the growth of new blood vessels to grow toward and eventually surround and penetrate the solid tumour. Angiogenesis is therefore a crucial component of solid tumour growth. Additionally, cancer cells actively invade the surrounding tissue through the over-expression and secretion of proteolytic enzymes such as matrix metalloproteases and urokinase plasminogen activators. These enzymes degrade the tissue locally and enable the cancer cells to both proliferate and actively migrate into the modified matrix. In this article, we present a new multiscale mathematical model for solid tumour growth which nonlinearly couples an improved model of tumor invasion with a model of tumour-induced angiogenesis. We perform simulations of the multi-scale model that demonstrate the importance of the nonlinear coupling between the development and remodelling of the vascular network, the blood flow through the network and the tumour progression. Consistent with clinical observations, the hydrostatic stress generated by tumor cell proliferation shuts down large portions of the vascular network dramatically affecting the flow, the subsequent network remodelling, the delivery of nutrients to the tumour and the subsequent tumor progression. In addition, extracellular matrix degradation by tumour cells is seen to have a dramatic affect on both the development of the vascular network and the growth response of the tumour. In particular, the newly developing vessels tend to encapsulate, rather than penetrate, the tumour and are thus less effective in delivering nutrients. Keywords solid tumour, avascular growth, angiogenesis, vascular growth, multi-scale mathematical model 1. Introduction Cancer growth, and as a particular example in this paper, solid tumour growth, is a complicated phenomenon involving many inter-related processes across a wide range of spatial and temporal scales, and as such presents the mathematical modeller with a correspondingly complex set of problems to solve. The aim of this paper is to formulate a multi-scale mathematical model of solid tumour growth, incorporating three key features: the avascular growth phase, the recruitment of new blood vessels by the tumour (angiogenesis) and the vascular growth and host tissue invasion phase. Solid tumours are known to progress through two distinct phases of growth the avascular phase and the vascular phase. The initial avascular growth phase can be studied in the laboratory by culturing cancer cells in the form of three-dimensional multicell spheroids. It is well known that these spheroids, whether grown from established tumour cell lines or actual in vivo tumour specimens, possess growth kinetics which are very similar to in vivo solid tumours. Typically, these avascular nodules grow to a few millimetres in diameter. Cells towards the centre, being deprived of vital nutrients, die and give rise to a necrotic core. Proliferating cells can be found in the outer cell layers. Lying between these two regions is a layer of quiescent (or hypoxic) cells, a proportion of which can be recruited into the outer layer of proliferating cells. Much experimental data has been gathered on the internal architecture of spheroids, and studies regarding the distribution of vital nutrients (e.g. oxygen) and metabolites within the spheroids have been carried out. See, for example, the recent reviews by Walles et al. (2007), Kim (2005), Kunz-Schughart et al. (2004), Chomyak and Sidorenko (2001) and the references therein. The transition from the relatively harmless and confined dormant avascular state to the vascular state, wherein the tumour possesses the ability to invade surrounding tissue and metastasise to distant parts of the body, depends upon the ability of the tumour to induce new blood vessels from the surrounding tissue to sprout towards and then gradually surround and penetrate the tumour, thus providing it with an adequate blood supply and microcirculation. Tumour-induced angiogenesis, the process by which new blood vessels develop from an existing vasculature, through endothelial cell sprouting, proliferation and fusion, is therefore a crucial part of solid tumour growth. Sustained angiogenesis is a hallmark of cancer (Hanahan and Weinberg 2000). Mature endothelial cells are normally quiescent and, apart from certain developmental processes (e.g. embryogenesis and wound healing), angiogenesis is generally a pathological process implicated in arthritis, some eye diseases and solid tumour development, invasion and metastasis. Tumour-induced angiogenesis is believed to start when a small avascular tumour exceeds a critical diameter (~2 mm), above which normal tissue vasculature is no longer able to support its growth. At this stage, the tumour cells lacking nutrients and oxygen become hypoxic. In response, the tumour cells secrete a number of diffusible chemical substances tumour angiogenic factors (TAF) - into the surrounding tissues and extracellular matrix (ECM). The TAF diffuses into the surrounding tissue and eventually reach the endothelial cells (EC) that line nearby blood vessels. ECs subsequently respond to the TAF concentration gradient by degrading the basement membrane surrounding the parent vessel, forming sprouts, proliferating and migrating towards the tumour. It takes approximately 10 to 21 days for the growing network to link the tumour to the parent vessel, and this vascular connection subsequently provides all the nutrients and oxygen required for continued tumour growth. An excellent summary of all the key cell-biological processes involved in angiogenesis can be found in the comprehensive review article of Paweletz and Knierim (1989). See also the recent review by Carmeliet (2005). Once vascularized the solid tumours grow rapidly as exophytic masses. In certain types of cancer, e.g. carcinoma arising within an organ, this process typically consists of columns of cells projecting from the central mass of cells and extending into the surrounding tissue area. The local spread of these carcinoma often assume an irregular jagged shape. By the time a tumour has grown to a size whereby it can be detected by clinical means, there is a strong likelihood that it has already reached the vascular growth phase. Cancers also possess the ability to actively invade the local tissue and then spread throughout the body. Invasion and metastasis are the most insidious and life-threatening aspects of cancer (Liotta and Stetler-Stevenson, 1991; Liotta and Clair, 2000). Indeed, the prognosis of a cancer is primarily dependent on its ability to invade and metastasize. Many steps that occur during tumour invasion and metastasis require the regulated turnover of extracellular matrix (ECM) macromolecules. The breakdown of these barriers is catalyzed by proteolytic enzymes released from the invading tumour. Most of these proteases belong to one of two general classes: many are metalloproteases (Parsons et al., 1997), while others are serine proteases (Andreasen et al., 1997; 2000). Proteases give cancers their defining characteristic - the ability of malignant cells to break out of tissue compartments. Motility, coupled with regulated, intermittent adhesion to the extracellular matrix and degradation of matrix molecules, allows an invading cell to move through the extracellular matrix (Liotta and Stetler-Stevenson, 1991; Lauffenburger and Horwitz, 1996; Friedl and Wolf 2003). However, proteolytic degradation of the extracellular matrix is essential for the key processes involved in tissue remodelling as well. These processes take place in a number of distinct physiological events in the healthy organism, such as trophoblast invasion, mammary gland involution, angiogenesis and wound healing. The most significant turning point in cancer, however, is the establishment of metastasis. The metastatic spread of tumour cells is the predominant cause of cancer deaths, and with few exceptions, all cancers can metastasize. Metastasis is defined as the formation of secondary tumour foci at a site discontinuous from the primary tumour (Liotta and Stetler-Stevenson, 1991; Liotta and Clair, 2000). Metastasis unequivocally signifies that a tumour is malignant and this is in fact what makes cancer so lethal. In principal, metastases can form following invasion and penetration into adjacent tissues followed by dissemination of cells in the blood vascular system (hematogeneous metastasis) and lymphatics (lymphatic metastases). Metastases can appear shortly after surgery but can also remain undetected for more than a decade before manifesting themselves clinically (King, 2000; Chambers et al., 2002; Fidler, 2002). This indicates that disseminated cancer cells can persist in a dormant state (either individually or as an avascular tumour spheroid), unable to form a progressively increasing tumor mass (Chambers et al., 2002). Such heterogeneity of outcome indicates that the fate of tumour cells that disseminate to distant organs before surgery must be regulated by either inherent cancer cell properties or the milieu of the target organs, or both. Identifying the mechanisms that keep metastases in their dormant, occult state is one of the most challenging and important avenues of cancer research (Chambers et al., 2002; Fidler, 2002). Since the seminal work of Greenspan (1976) the mathematical modelling of avascular solid tumour growth, like its subject, has been rapidly expanding. Most models in this area consist of systems of nonlinear partial differential equations, and may be described as macroscopic. The review paper of Araujo and McElwain (2004) provides an excellent overview. See also the recent reviews by Quaranta et al. (2005), Byrne et al. (2006), Sanga et al. (2006), Graziano and Preziosi (2007) and Roose et al. (2007). Likewise, modelling tumour-induced angiogenesis has a well-established history beginning with the work of Balding and McElwain (1983). The review paper of Mantzaris et al. (2004) again providing an excellent overview of the work in this area. However, unlike avascular growth and angiogenesis, vascular tumour growth has received considerably less attention in the mathematical modelling literature. Recently, Zheng et al. (2005) developed and coupled a level-set method for solid tumour growth with a hybrid continuous-discrete model of angiogenesis originally developed by Anderson and Chaplain (1998). This work served as a building block for studies of chemotherapy (Sinek et al. (2004)) and morphological instability and tumour invasion (Cristini et al. (2005); Frieboes et al. (2006)). Hogea et al. (2006) have also begun to investigate tumour induced angiogenesis and vascular growth using a level-set method coupled with a continuous model of angiogenesis. Following the strategy pioneered by Zheng et al., Frieboes et al. (2007) coupled a mixture model with a lattice-free continuous-discrete model of angiogenesis (originally developed by Planck and Sleeman (2004)) and studied vascular tumour growth in three dimensions. In these works, however, the effects of blood flow through and subsequent remodelling of the vascular network were not included. Recently, the effects of blood flow through a vascular network on tumour growth were considered by Alarcon et al. (2005), Lee et al (2006), Bartha and Rieger (2007), Welter et al. (2007) using cellular automaton (CA) tumour growth models coupled with network models for the vasculature. These authors investigated vascular network inhomogeneities, the stress-induced collapse of blood vessels and the implications for therapy. Because of the computational cost of simulating cell growth using CA, these studies are limited to small scales. In this paper we couple an improved continuum model of solid tumour invasion (following Macklin and Lowengrub 2007) that is capable of spanning the 102mm-cm scale and accounts for cell-cell, cell-ECM adhesion, ECM degradation, and tumor cell migration, proliferation, and necrosis with a model of tumour-induced angiogenesis (following McDougall et al. (2006)) that accounts for blood flow through the vascular network, non-Newtonian effects and vascular network remodelling, due to wall shear stress and mechanical stresses generated by the growing tumour, to produce a new multi-scale model of vascular solid tumour growth. As in Zheng et al. (2005), the invasion and angiogenesis models are coupled through the tumor angiogenic factors (TAF) that are released by the tumor cells and through the nutrient extravasated from the neo-vascular network. As the blood flows through the neo-vascular network, nutrients (e.g. oxygen) are extravasated and diffuse through the ECM triggering further growth of the tumour, which in turn influences the TAF expression. In addition, the extravasation is mediated by the hydrostatic stress generated by the growing tumour and, as mentioned above, the hydrostatic stress also affects vascular remodelling by restricting the radii of the vessels. The vascular network and tumour progression are also coupled via the ECM as both the tumour cells and ECs upregulate matrix degrading proteolytic enzymes which cause localized degradation of the ECM which in turn affects haptotactic migration. We perform simulations of the multi-scale model that demonstrate the importance, on tumour invasion of the host tissue, of the nonlinear coupling between the growth and remodelling of the vascular network, the blood flow through the network and the tumour progression. Consistent with clinical observations, the hydrostatic stress generated by tumour cell proliferation shuts down large portions of the vascular network dramatically affecting the flow, the subsequent network remodelling, the delivery of nutrients to the tumour and the subsequent tumour progression. In addition, ECM degradation by tumour cells is seen to have a dramatic affect on both the development of the vascular network and the growth response of the tumour. In particular, when the ECM degradation is significant, the newly formed vessels tend to encapsulate, rather than penetrate, the tumour and are thus less effective in delivering nutrients. The outline of the paper is as follows. In section 2, the mathematical models are presented. In section 3, the numerical methods are briefly described. In section 4, numerical results are presented and conclusions and future work are discussed in section 5. Details of the mathematical modelling and numerical methods are presented in the appendices. 2. The Mathematical Model In this exposition, we present the non-dimensional model, starting first with the model of tumour invasion in section 2.1 and followed by the model of tumour-induced angiogenesis in section 2.2. Here, time is non-dimensionalized by the characteristic tumor cell proliferation time (i.e., EMBED Equation.DSMT4 where  EMBED Equation.DSMT4 is the mitosis rate) and space is non-dimensionalized by the characteristic diffusion penetration length (i.e.,  EMBED Equation.DSMT4 , where EMBED Equation.DSMT4 and  EMBED Equation.DSMT4 are characteristic values of the oxygen diffusion coefficient and uptake rate in the proliferating tumour region respectively). In the following, quantities defined with an overbar correspond to non-dimensional constants. The non-dimensionalization of the parameters and the corresponding values used in the numerical simulations are presented in Tables 4-6. For the non-dimensionalization and further biological discussion of the models, the reader is referred to Macklin and Lowengrub (2006, 2007) and McDougall et al. (2006). 2.1 The tumour invasion model To accurately model tumour growth in heterogeneous tissues, we develop a mathematical model that accounts for spatially dependent cell-necrosis, cell-apoptosis, cell-cell and cell-matrix adhesion, matrix degradation, cell-proliferation and cell-migration. The model is based on continuum reaction-diffusion equations that describe these processes and is a generalization and improvement of earlier models (see the reviews listed previously and recent work by Macklin and Lowengrub 2007b, 2007c). We present the model in two-dimensions but it is equally valid for the three-dimensional case as well. Let  denote a tumour mass, and let  denote its boundary. The tumour can be divided into three regions: a proliferating rim P where the tumour cells have sufficient nutrient levels for proliferation; a hypoxic/quiescent region H where the nutrient levels are too low for normal metabolic activity but not so low that the cells begin to die; and a necrotic region N where the nutrient level has dropped so low that the tumour cells have died and are degraded. Because necrosis is irreversible, we track the necrotic core and its interface N separately of the tumour interface. See Figure 1. 2.1.1 Nutrient Transport We model the net effect of nutrients (e.g., oxygen and glucose) and growth-promoting and  inhibiting factors with a single nutrient . Here, we focus our attention on the role of oxygen which is supplied by the vascular network via the red blood cells. This can be modelled using the haematocrit which represents the volume fraction of red blood cells contained in the blood. Oxygen, and other nutrients, are supplied by the pre-existing bulk vasculature and the neo-vasculature at rates bulk and neo, diffuses throughout the cancerous and non-cancerous tissue, is uptaken in the non-necrotic portions of the tumour, and decays elsewhere (see below). Wherever the oxygen level inside the tumour drops below a threshold value H, the tumour cells become hypoxic (quiescent), cease proliferating and uptake nutrient at a lower rate. If the oxygen level falls further below a threshold value N, then the tumour cells become necrotic. Inside the necrotic core, oxygen reacts with cellular debris to form reactive oxygen species (Kloner and Jennings 2001, Galaris et al. 2006), which we model by a decay term. Since oxygen diffusion occurs more rapidly than cell-mitosis (the time scale on which the equations are non-dimensionalized), these processes are described by the quasi-steady reaction-diffusion equation  EMBED Equation.3 , (1) where D is the diffusion coefficient, the parameter  combines the effects of oxygen uptake and decay and takes the form  EMBED Equation.3 , (2) where ps and q are smooth interpolating functions (the precise forms are given in Appendix 3) and E0 the density of the original ECM which is used to assess changes in uptake/decay in the host microenvironment (see section 2.1.3 below). The interpolating function  EMBED Equation.DSMT4 , where sH and sN are the oxygen concentration thresholds for quiescence and necrosis respectively and  EMBED Equation.3  is the rate of oxygen uptake by quiescent cells in the hypoxic tumour. Further,  EMBED Equation.3  and  EMBED Equation.3  are the rates of oxygen uptake in the host microenvironment and in the proliferating tumour regions respectively, and  EMBED Equation.3  is the rate of oxygen decay in the necrotic portion of the tumour. We note that because the location of the viable, hypoxic, and necrotic tumour regions depends upon the past and present values of the oxygen level , the uptake/decay term  introduces nonlinearity. The two remaining sources  EMBED Equation.DSMT4  and  EMBED Equation.DSMT4  in Eq. (1) reflect the oxygen-tissue transfer from the pre-existing and neo-vascular blood vessels respectively and are given by:  EMBED Equation.3  (3) and  EMBED Equation.3 , (4) where  EMBED Equation.DSMT4  and  EMBED Equation.DSMT4 are constant transfer rates from the pre-existing and neo- vessels. Here, Bpre is the (non-dimensional) blood vessel density of the pre-existing vessels whose locations are assumed to be unchanging in time. In fact, we take a uniform distribution of pre-existing vessels in the host tissue and Bpre satisfies Eq. (19a) below where MDE is assumed to degrade the pre-existing vasculature. The function EMBED Equation.DSMT4 is the characteristic or indicator function of the neo-vessels (i.e., equal to 1 at the locations of the neo-vessels), and 1W is the characteristic function of the tumour region W (i.e., equal to 1 inside the tumor and is 0 in the tumor exterior). Further, P is the oncotic (solid/mechanical) pressure, Pvessel and h are the dimensional pressure and the haematocrit in the neo-vascular network, respectively. The constants  EMBED Equation.DSMT4 and  EMBED Equation.DSMT4 reflect the normal value of haematocrit in the blood (generally about 0.45) and the minimum haematocrit needed to extravasate oxygen, respectively. The haematocrit is modelled via the blood flow in the vascular network and is determined from the angiogenesis model. This provides one aspect of the coupling between the tumour growth and angiogenesis models. A second mode of coupling between the two models occurs through the cutoff function c(Pvessel,P) which is given by:  EMBED Equation.DSMT4 , (5) where pcutoff is a cubic, interpolating polynomial given in Appendix 3. Namely, large oncotic pressures may prevent extravasation and transfer of oxygen from the vessels into the tissue. Later, we will discuss how the oncotic pressure may also constrict the neo-vessels. Further, in Eq. (5),  EMBED Equation.DSMT4  , (6) where EMBED Equation.DSMT4  is a characteristic pressure scale and  EMBED Equation.DSMT4 is a scale factor. Note that we could have analogously taken the oxygen transfer rate from the pre-existing vessels to also be coupled to the haematocrit and blood vessel pressure. This will be explored in a future work. The oxygen source terms in E qs. (3) and (4) are designed such that for sufficiently large transfer rates EMBED Equation.DSMT4  and  EMBED Equation.DSMT4 , the oxygen concentration  EMBED Equation.DSMT4 at the spatial locations of the pre-existing and neo- vessels. In practice, we will take  EMBED Equation.DSMT4  large but  EMBED Equation.DSMT4 small which models the supply of only a small amount of oxygen in the host tissue from pre-existing vessels. We will assume a parent vessel, located at the boundary of the computational microenvironment domain as discussed below, supplies the bulk of the oxygen in the host tissue. Note that oxygen flux conditions across the pre-existing and neo-vessels could be imposed (e.g., Alarcon et al. 2005). The boundary conditions for Eq. (1) are taken to be a combination of Dirichlet, Neumann conditions. In particular, in the simulations we present below, we assume a parent vessel coincides the upper boundary of the computational domain. Therefore, a Dirichlet condition, s =1, is posed along the upper boundary. Zero Neumann conditions, EMBED Equation.DSMT4 , are imposed along the other boundaries of the computational domain. 2.1.2 Tumour mechanics and the cell-velocity The tumor cells, the ECM and host, non-cancerous cells are influenced by a combination of forces which contribute to the cellular velocity field. The proliferating cells generate an oncotic mechanical pressure (hydrostatic stress) that also exerts force on the ECM and host cells. The cells respond to pressure variations by overcoming cell-cell and cell-ECM adhesion and migrating through the microenvironment. The ECM may also deform, degrade and remodel in response to pressure and to the chemical factors released by the cells. The cells may respond haptotactically to adhesion gradients in the ECM. Following previous work, we assume that all solid phases move with a single cellular velocity field and we model the cellular motion within the ECM as incompressible fluid flow in a porous medium. In the future, we plan to use mixture models (e.g., Ambrosi and Preziosi 2002, Byrne and Preziosi 2003, Araujo and McElwain 2005a, 2005b) to relax these assumptions. In this simplified description of tumour mechanics used here, Darcys law is taken as the constitutive assumption and thus the velocity is proportional to the forces in the problem. See Ambrosi and Preziosi (2002) and Byrne and Preziosi (2003) for a motivation of this approach from a mixture modelling perspective. Accordingly, the non-dimensional velocity is given by  EMBED Equation.3 , (7) where EMBED Equation.DSMT4  is the cell-mobility which models the net effects of cell-cell and cell-matrix adhesion, E is the ECM density (e.g. a non-diffusible matrix macromolecule such as fibronectin, collagen or laminin) and EMBED Equation.DSMT4  is the haptotaxis coefficient. Models for m and cE are given below in section 2.1.2. Further assuming that the density of tumour cells is constant in the viable region, the growth of the tumour is then associated with the rate of volume change:  EMBED Equation.3 , (8) where EMBED Equation.DSMT4  is the non-dimensional net proliferation rate. This implies that the non-dimensional pressure satisfies:  EMBED Equation.3 . (9) We assume that in the proliferating region, cell-mitosis is proportional to the amount of nutrient present and that apoptosis may occur. Volume loss may occur in the necrotic core and there is no proliferation in either the host microenvironment or the hypoxic/quiescent regions. We therefore take  EMBED Equation.3 , (10) where A is the non-dimensional apoptosis rate (pre-programmed cell death); and GN is the non-dimensional rate of volume loss in the necrotic core as water is removed and cellular debris are degraded. Assuming a uniform cell-cell adhesion throughout the tumour, cell-cell adhesion can be incorporated as a surface-tension like jump boundary condition at the tumour-host interface  EMBED Equation.DSMT4 :  EMBED Equation.3 , (11) where G is a non-dimensional parameter that measures the aggressiveness of the tumour (the strength of cell proliferation relative to cell-cell adhesion) and  EMBED Equation.DSMT4 is the mean curvature of the interface. At the necrotic boundary SN we assume P is continuous. We assume that no voids form and therefore we take  EMBED Equation.DSMT4  which implies that  EMBED Equation.DSMT4 , (12) where n is the unit outward normal to S. For simplicity, we will also assume that EMBED Equation.DSMT4 . At necrotic boundary, we assume analogous conditions. The velocity of the tumour-host interface EMBED Equation.DSMT4  is then given by:  EMBED Equation.3 , (13) and the velocity of the necrotic boundary  EMBED Equation.DSMT4  is  EMBED Equation.3 , (14) where nN is the outward unit normal vector along N. In the far-field at the boundaries of the computational domain, the pressure is assumed to satisfy zero Neumann boundary conditions  EMBED Equation.DSMT4 . 2.1.3 Tumour-Microenvironment Interaction We model tumour microenvironment by introducing an extracellular matrix density E that represents the density of non-diffusible matrix macromolecules such as fibronectin, collagen, elastin and laminin etc. In addition, as mentioned earlier, we keep track of the density E0 of the original ECM and the pre-existing blood vessel density Bpre to assess the level of oxygen uptake and supply, respectively, in the microenvironment. The tumour interacts with the microenvironment by responding to the nutrients supplied by the pre-existing and the neo-vasculature (e.g. see Eq. (1)), remodeling the ECM locally by secreting both MDE and ECM macromolecules and by a hetereogeneous response to pressure and ECM adhesion gradients through non-constant cell-mobility and haptotaxis coefficients. In order for tumours cells to migrate into the porous matrix, they must overcome cell-matrix adhesion. However, in experiments, a maximum migration speed is obtained that depends on the level of integrin expression (e.g. Pelacek et al. (1997), DiMilla et al. (1991)) and correspondingly a non-monotonic dependence of cell migration velocity on integrin expression and adhesion gradients in the ECM has been predicted (DiMilla et al, Dickinson and Tranquillo (1993)). This has been explained by the fact that while some integrins are required for focal adhesion based migration, too much focal contact strength can retard the detachment of cells trailing edge from the ECM. While we do not model integrin expression directly here, we take this effect into account by making the haptotaxis coefficient a non-monotone function of E:  EMBED Equation.3 , (15) where  EMBED Equation.DSMT4  are the non-dimensional haptotaxis in low/high-density ECM, pc is a non-monotone interpolating function with a maximum  EMBED Equation.DSMT4 located at  EMBED Equation.DSMT4 . See Appendix 3 for the precise form of pc. Although the mobility m may also be non-monotone, for simplicity, we take a monotone decreasing function of E here:  EMBED Equation.3 , (16) where pm is a smooth interpolating function (see Appendix 3 for a precise form). In a future work, we will investigate non-monotonic cell mobilities m. In addition, the mobility and chemotaxis parameters may also be functions of oxygen concentration s as hypoxic conditions may result in upregulation of HIF-1 alpha target genes that may result in decreased cell-cell and cell-matrix adhesion, among other effects, and therefore enable cells to more easily migrate through and invade the tumour microenvironment (e.g., see Kaur et al (2005), Erler et al. (2006), Pouyssegur et al. (2006)). These effects will also be explored in a forthcoming work. In order to migrate through the ECM and invade the host tissue, tumour cells secrete matrix degrading proteolytic enzymes (MDE), e.g. matrix metalloproteases and urokinase plasminogen activators, which cause the degradation of the ECM, provide space for the cells, and enhance the attachment of the cells to ECM macromolecules enabling the cells to exert traction forces to propel themselves through the ECM. In addition, the tumour cells remodel the ECM by secreting insoluble matrix macromolecules and possibly re-orienting them. We note that during the angiogenetic response of the host vasculature, an analogous molecular cascade occurs as tumour angiogenesis factors (TAF) and ECM macromolecules (e.g. fibronectin, collagen, laminin) bind to specific membrane receptors on ECs and activate the cells migratory machinery. This leads to a remodelling of ECM similar to that described above for tumour cells. Here, we will not consider the effect of orientational re-ordering. We model the remaining processes as follows. For the MDE, we take  EMBED Equation.3  (17) where M is the nondimensional MDE concentration,  EMBED Equation.DSMT4  is the diffusion coefficient (assumed to be constant),  EMBED Equation.DSMT4  and EMBED Equation.DSMT4 are the non-dimensional rates of production of MDE by the viable tumour cells ( EMBED Equation.DSMT4 ) and the sprout-tip ECs respectively. Further,  EMBED Equation.DSMT4 is the rate of decay (it is assumed that MDE is not used up as a result of the interaction with the ECM (Quaranta, private communication)). Finally,  EMBED Equation.DSMT4 is the characteristic function of the sprout tips. Because the diffusion coefficient of MDE, DM, is much smaller than that for oxygen diffusion the full time-dependent diffusion equation is used (Stephanou et al., 2006). In the far-field (boundary of the computational domain), we take the zero Neumann boundary conditions EMBED Equation.3 . The ECM density satisfies:  EMBED Equation.3 , (18) where EMBED Equation.DSMT4 and EMBED Equation.DSMT4  are the non-dimensional rates of production of ECM by the viable tumour cells and sprout-tip ECs and the EMBED Equation.DSMT4 is the non-dimensional rate of matrix degradation by the MDE. Finally, the original ECM and the pre-existing blood vessel density are assumed to be degraded by the MDE:  EMBED Equation.DSMT4 and  EMBED Equation.DSMT4  , (19) where  EMBED Equation.DSMT4 and  EMBED Equation.DSMT4 are non-dimensional degradation rates. 2.1.4 Tumour Angiogenic Factors When tumour cells become hypoxic/quiescent, they are assumed to secrete tumour angiogenic factors (TAF), which diffuse into the surrounding tissue and attract ECs. ECs respond to the TAF by binding with it, proliferating and chemotaxing up the TAF gradient. The diffusion coefficient of TAF is similar to that of oxygen and so we model the production, diffusion, decay, and binding of TAF by  EMBED Equation.3  (20) where T is the non-dimensional TAF concentration,  EMBED Equation.DSMT4 is the diffusion coefficient (assumed to be constant) and  EMBED Equation.DSMT4 ,  EMBED Equation.DSMT4 , and  EMBED Equation.DSMT4 denote the non-dimensional production, natural decay and binding rates of TAF. In the far-field at the boundary of the computational domain, we also take zero Neumann boundary conditions EMBED Equation.3 . 2.2 Angiogenesis model We begin with a description, in section 2.2.1, of an initial mathematical model for the growth of a hollow capillary network in the absence of any blood flow. This follows Anderson and Chaplain (1998). Then, following McDougall et al. (2006), we will add the effects of blood flow and vascular network remodelling in sections 2.2.2 and 2.2.3, respectively. 2.2.1 Basic network model As described earlier, TAF and ECM macromolecules bind to specific membrane receptors on ECs and activate the cells migratory machinery. The model of EC migration given below describes how capillary sprouts emerging from a parent vessel migrate towards a tumour, leading to the formation of a vascular network that supplies nutrients for continued development. See figure 2. At this level, since there is no flow or vessel remodelling, this model may perhaps be considered more appropriate at describing in vitro endothelial cell migration and capillary sprout formation. The model, inspired by the tumor angiogenesis model developed by Anderson and Chaplain (1998) assumes that endothelial cells migrate through (i) random motility, (ii) chemotaxis in response to TAF released by the tumour and (iii) haptotaxis in response to ECM gradients. If we denote by n the non-dimensional endothelial cell density per unit area, then the non-dimensional equation describing EC conservation is given by  EMBED Equation.DSMT4 . (21) See McDougall et al. (2006) and Appendix 3 for the non-dimensionalization. The diffusion (random migration) coefficient is  EMBED Equation.DSMT4  (assumed to be constant), and the chemotactic and haptotactic migration are characterised by the functions EMBED Equation.DSMT4 , which reflects the decrease in chemotactic sensitivity with increased TAF concentration and EMBED Equation.DSMT4 , where for simplicity we have taken the haptotactic migration parameter to be constant. In a future work, we will investigate the heterogeneous response of the ECs to the ECM and oxygen gradients as discussed earlier in section 2.1.3. The coefficients EMBED Equation.DSMT4 ,  EMBED Equation.DSMT4  and  EMBED Equation.DSMT4 characterise the non-dimensional random, chemotactic and haptotactic cell migration respectively. The displacement of each individual EC, located at the tips of each growing sprout, is given by the discretised form of the EC mass conservation equation (21) on a regular Cartesian mesh. The migration of each cell is consequently determined by a set of coefficients (P0-P4) emerging from this equation, which relate to the likelihood of the cell remaining stationary, moving left, right, up or down. These coefficients incorporate the effects of random, chemotactic and haptotactic movement and depend upon the local chemical environment (ECM density and TAF concentration). Proliferation of the endothelial cells at the capillary tips and branching at capillary tips are implemented in the model at the discrete level. Tip branching depends on the TAF concentration at a given spatial location (see Table 1 below and Anderson and Chaplain (1998) for details.). Using the above model it is possible to generate hollow capillary networks which are structurally similar to those observed experimentally. 2.2.2 Modelling blood flow in the developing capillary network Blood is a complex multiphase medium, composed of many different constituents, including: red blood cells (erythrocytes), white blood cells (leukocytes), and platelets involved in clotting cascades. These solid elements represent approximately 45% of the total blood composition red cells are predominant and are carried in the plasma, which constitutes the fluid phase. A measure of the solid phase is given by the blood haematocrit, which represents the volume fraction of red blood cells contained in the blood. The average human haematocrit has a value of around 0.45. Because of its multiphasic nature, blood does not behave as a continuum and the viscosity measured while flowing at different rates in microvessels is not constant. The direct measurement of blood viscosity in living microvessels is very difficult to achieve with any degree of accuracy. However, by comparing the flow distribution in a numerical network (generated by a mathematical model) with a similar experimental system, Pries et al. (1996) determined a relationship between the apparent viscosity of blood, the blood haematocrit and the radius of the vessel, though which the blood is flowing, that provides a good fit of microvascular experimental data. The relationship is given as follows:  EMBED Equation.DSMT4  (22) where mplasma is the plasma viscosity, mrel is the relative viscosity, EMBED Equation.DSMT4  is the (nondimensional) viscosity corresponding to the normal value of the discharge haematocrit (i.e.,  EMBED Equation.DSMT4 ), R is the dimensional vessel radius (in mm), R* is a radius scale factor (taken to be equal to 1 mm), and f(h,R) is a modulating function of the haematocrit and vessel radius. These effects are modelled by:  EMBED Equation.DSMT4  (23) The apparent blood viscosity (e.g. Eq. (22a)) generally increases with decreasing capillary radius, although the precise relationship is nonlinear since it is actually haematocrit-dependent. In order to calculate flow within the entire interconnected network of capillaries, it is first necessary to decide upon a local relationship between the pressure gradient  EMBED Equation.DSMT4  and flow rate  EMBED Equation.DSMT4  at the scale of a single capillary element of length L and radius R. Such a relationship in the case of a non-Newtonian fluid can be approximated by the following Poiseuille-like expression:  EMBED Equation.DSMT4 . (24) In order to determine the pressure (and flow rate) and in the vascular network of interconnected capillary elements having distributed radii, one simply conserves mass (or flow if the fluid is incompressible) at each junction where capillary elements meet. See Figure 3. Hence, for each node the following expression can be written:  EMBED Equation.DSMT4  (25) where the index k refers to adjacent nodes and N=4 in a fully connected regular 2D grid as considered in this paper (or N=6 in 3D). This procedure leads to a set of linear equations for the nodal pressures (Pvessel,i) which can be solved numerically using any of a number of different algorithms including successive over-relaxation (SOR). Once the nodal pressures are known, Eq. (24) can be used to calculate the flow in each capillary element in turn. A more complete discussion of the procedure can be found in McDougall et al. (2002). The evolution of haematocrit h in the vessels is also calculated using mass conservation once the flow is determined. 2.2.3. Capillary vessel adaptation and remodelling Blood rheological properties and microvascular network remodelling are interrelated issues, as blood flow creates stresses on the vascular wall (shear stress, pressure, tensile stress) which lead to adaptation of the vascular diameters via either vasodilatation or constriction. In turn, blood rheology (viscosity, haematocrit, etc.) is affected by the new network architecture consequently, we should expect adaptive angiogenesis to be a highly dynamic process. We follow the work of Pries et al. (1995, 1996, 1998, 2001a) in incorporating vessel adaptation into our model. In particular, we consider a number of stimuli that affect the vessel diameters. In particular, we account for the influence of the wall shear stress (Swss), the intravascular pressure (Sp), a metabolic mechanism depending on the blood haematocrit (Sm), as well as the natural tendency for vessels to shrink (Ss). These stimuli form a basic set of requirements in order to obtain stable network structures with realistic distributions of vessels diameters and flow velocities. The theoretical model for vessel adaptation assumes that the change in a flowing vessel radius DR over a time step Dt, where time is scaled by the rate of the response of the vessel to wall shear stress ( EMBED Equation.DSMT4 ), is proportional to both the global stimulus acting on the vessel and to the initial vessel radius R, i.e.  EMBED Microsoft Equation 3.0 , (26) We refer the reader to Appendix 1 for the definitions of the stimuli and a brief discussion. More details may be found in McDougall et al. (2006). After the radius of the vessel is updated according to Eq. (26), the effect of the oncotic mechanical pressure P, generated by the proliferating and invading tumour, on the vessel radius is then taken into account. The tendency of the oncotic pressure to shrink the vessel is modelled by the simple cutoff:  EMBED Equation.DSMT4 , (27) where c(P,Pvessel) is the cutoff function introduced earlier in Eq. (5) and Rmin is a threshold minimum radius. This provides another means of coupling tumour invasion (and mechanics) with the angiogenic response and the developing neo-vascular network. In particular, the solid/mechanical pressure may constrict and cut off vessels in the neovasculature. To prevent singularities in practice, the radius of the vessel is constrained to lie between 2.0mm and 14mm which is the size of the parent capillary. Inclusion of the above mechanisms into our modelling framework now allows us to simulate dynamic remodelling of a flowing vasculature. This significant improvement in angiogenesis modelling, introduced by McDougall et al. (2006), allows us to describe vascular growth in a far more realistic manner, with areas of the capillary network dilating and constricting in response to variations in perfusion-related stresses, stimuli and pressure mechanical forces exerted on the host microenvironment by the invading tumour. The final step in the development of the complete dynamic adaptive tumour-induced angiogenesis (DATIA) model is to couple the network flow modelling approach outlined in this section to the hollow capillary model derived from the endothelial cell migration equations described earlier. This is achieved through the role of wall shear stress. Wall shear stress is known to play a leading role in the growth and branching of capillary vessel networks Pries et al. (2001a, 2001b). In order to bring the morphological and the physiological concepts together (Thompson, 1917), the cell migration and flow models are coupled by incorporating the mechanism of shear-dependent vessel branching in addition to sprout-tip branching via local TAF concentrations. This enables the capillary network structures to adapt dynamically through adjuvant vessel branching in areas of the network experiencing increased shear stresses following anastomosis elsewhere in the system. We note that because the shear stress is due to the blood flowing through the capillaries, vessel branching can only occur after some degree of anastomosis has taken place. Therefore, the early stages of angiogenesis are primarily characterised by branching at the capillary tips which depends only on the TAF concentration. The combined effects of the local wall shear stress and TAF concentration upon vessel branching probability have been implemented in the model as described in Table 1. In the absence of quantitative experimental data, the probabilities chosen for the vessel branching process have been defined on a qualitative basis and reflect the combined influence of the wall shear stress (WSS) and local TAF concentration. High values of WSS in tandem with high local TAF concentrations lead to a higher branching probability, whilst lower values of one or both of WSS and TAF concentration lead to lower branching probability. For each range of WSS (linearly distributed in [0,1]), the corresponding TAF probability profile has been obtained via a linear scaling of the values reported in McDougall et al. (2002) and Stphanou et al. (2005a, 2005b). As mentioned above, in the absence of WSS, TAF-dependent sprout tip branching is the only means by which a migrating vessel can bifurcate. Sprout tip branching is performed using the algorithm developed by Anderson and Chaplain (1998) and the corresponding tip branching probabilities are shown in Table 2. 3. Numerical Schemes 3.1 Tumour invasion model The tumour invasion model described in Section 2.1 consists of a coupled system of nonlinear, elliptic and parabolic (reaction-diffusion) differential equations that must be solved on a complex, moving domain where the motion of the tumour/host boundary depends on gradients of the solutions to these equations. Further, one of these solutionsthe pressureis discontinuous across the tumour/host interface where the discontinuity depends on the geometry (i.e. the curvature) of the interface which is an additional source of nonlinearity. Therefore, standard finite difference methods cannot be used to accurately solve the system. Instead, specialized methods that can accurately take into account discontinuities in solutions and complex domains must be used. Here, we use a ghost-cell/level-set method and adapt and extend the numerical techniques we recently developed (Macklin and Lowengrub 2005, 2006, 2007) to solve this system. In this approach, the equations are discretized on a regular Cartesian mesh and the difference stencils near discontinuities are modified. We note that other alternatives exist (see the discussion in Macklin and Lowengrub 2007), but an advantage of our approach is that it can be implemented in a dimension-by-dimension manner, making the extension to 3D straightforward, and our algorithm is simpler to implement than the alternative approaches. In this approach, the interface is captured as the zero set of an auxiliary function (the level-set function)  satisfying  < 0 inside ,  > 0 outside , and  = 0 on the tumour/host interface . Typically  is taken to be an approximation to the signed distance function, i.e.  EMBED Equation.DSMT4 . See Figure 2. The interface normal and curvature can easily be calculated from . The interface N separating viable tumour cells from the necrotic cells is also captured using additional level set function boundary N that satisfies the same properties as , only with N and N in place of  and . Away from , the elliptic/parabolic equations can be discretized using centered finite differences. However, near the interface, the difference stencils need to be modified to account for possible jumps in solutions and in their normal derivatives. To do this, ghost cells on either side of the interface are introduced and the variables are extrapolated across the interface to ensure that the difference stencil effectively does not include nodes on the other side of the interface. The resulting nonlinear system is solved using an iterative algorithm. These techniques are described in Appendix 2. 3.2 The dynamic adaptive tumour-induced angiogenesis model For a fixed tumour geometry and TAF distribution, the tumour vasculature is first grown using the basic network model given in section 2.2.1; capillary tips may branch or anastomose during this stage. Further, the Cartesian mesh for the tumour growth system coincides with that used for the neo-vascular network. After a certain period of time, referred to as the capillary growth duration time, the fluid flow is solved in the fixed neo-vascular network and then the network is dynamically remodelled, following the algorithm described in sections 2.2.2 and 2.2.3 respectively. During the simulation of the flow, a CFL condition is imposed on the time step:  EMBED Equation.DSMT4 where  EMBED Equation.DSMT4  and  EMBED Equation.DSMT4 are the velocity and flow rate in a capillary element. The minimum is taken over the neo-vascular network. This ensures haematocrit remains conserved during the simulation (e.g., McDougall et al. 2002). Then process of blood flow, followed by remodelling, is repeated for an amount of time referred to as the flow duration time. 3.2 Overall Computational Solution Technique Initially, the avascular tumour, the pre-existing vascular network, the oxygen, ECM and MDE concentrations are given. We will consider a single parent vessel placed at the top of the computational domain (see Figs. 1 and 2). The algorithm then consists of iterating the following steps. (1). Solve Eq. (1) for the oxygen concentration where the oxygen source in Eq. (4) is obtained from the haematocrit and the pressure in the existing vascular network and the tumour mechanical pressure from the previous time step. We then use the solution s to update the position of the necrotic core:  EMBED Equation.3 , and to identify the hypoxic region H. As described above, the necrotic core is expanded to include previously necrotic tissue plus any tumour tissue where the oxygen level has dipped below the necrotic threshold N. We then rebuild N as a level set function that represents the updated region N. (Please see the Appendix 2, Macklin and Lowengrub (2007) and the level set references above for information on initializing a new level set function.) (2). Solve Eq. (20) for the tumor angiogenic growth factor (TAF) and update the MDE and ECM according to Eqs. (17) and (18) respectively. (3). Determine the cellular mobility and solve for the tumor biomechanical pressure from Eq. (9). (4). Update the position of the tumour/host interface  and the necrotic/viable N by advecting the level set functions  and N with the appropriate velocities as described in Appendix 2 (see Eq. (A2.3)). If necessary, the level-set functions are re-initialized via Eq. (A2.4) to be local distance functions to  and N. (5). From the updated tumour position, TAF, MDE and ECM fields, the neo-vascular network is grown using the basic network model. (6). The process (1)-(5) is repeated until the growth duration time interval is reached. At this point, the fluid flow in the neo-vascular network is determined and the network is adapted. The hydrostatic pressure P and the TAF are held fixed during this process. The flow and network adaption are repeated (for fixed tumour and capillary tip positions) until the flow duration time is reached. (7). Go to (1) and repeat the algorithm. 4. Computational Results In this work, we shall focus upon tumour growth coupled to angiogenesis in a square 4mm x 4mm region. Although we solve the non-dimensional equations, we present dimensional results using the length scale  EMBED Equation.DSMT4 and the time scale EMBED Equation.DSMT4 . A parent capillary vessel is located at the top of the computational domain. A pre-existing vasculature is assumed to exist and provides a small level of nutrient uniformly throughout the host tissue domain. Initially, a small cluster of proliferating cells is placed approximately 3mm from the parent vessel. The initial ECM is taken to be nearly constant (=1) but with small random perturbations uniformly distributed throughout the computational domain. See the time t=0 plot in Fig. 7(a). Accordingly, whenever we calculate gradients of E, we actually calculate the gradient of smoothed E where a Gaussian smoothing with standard deviation 3.0 is used (see Macklin and Lowengrub 2005). We begin by demonstrating that in the absence of tumour-induced angiogenesis, the small tumour cluster grows to an avascular tumour (2D) spheroid. Actually, since there is a pre-existing vasculature this is an abuse of notation, however, we still refer to this case as avascular since there is no neo-vascular network. Then, tumour-induced angiogenesis is initiated and we present several simulations of angiogenesis and vascular growth. Finally, we examine the effect of increased ECM degradation by MDE and its affect on avascular and vascular growth. The parameters, and non-dimensionalization, used in the simulations are given in Tables 4-6. 4.1 Avascular growth to a multicellular (2D) spheroid In Fig. 7(a), we present the growth of an avascular tumour. The spatial grid is  EMBED Equation.DSMT4  and the time step Dt=0.05 and is adapted to satisfy the Courant-Friedrichs-Lewy (CFL) condition (see Macklin and Lowengrub 2006, 2007). The red, blue and brown colors denote WP, WH, WN the proliferating, hypoxic/quiescent and necrotic regions respectively. The non-dimensional oxygen and ECM concentrations and the solid (oncotic) pressure are also shown. The oxygen diffuses only a short distance (about 0.2mm from the diffusion length) from the parent vessel as can be observed from the figure. However, the pre-existing vasculature (which yields a background oxygen concentration of approximately 0.4), provides enough oxygen for the tumour to grow. As the tumour grows, the pressure in the proliferating region increases, the oxygen is depleted in the tumour and the ECM is degraded. A hypoxic/quiescent core forms at about 9 days when the tumour radius is approximately 0.34mm. While the tumour continues to grow and degrade the extracellular matrix, the pressure decreases and the tumour growth starts to slow, as can be seen in Fig. 7(b). A necrotic core forms around day 15 when the radius of the tumour is approximately 0.5mm. The pressure drops significantly to reflect the volume loss in the necrotic core associated with the break-down of the necrotic cells and the growth of the tumour slows even further as the tumour approaches a steady state. As the growth of the tumour slows, the ECM degradation becomes more pronounced. This actually causes a competition between two effects: the pressure-induced motion, which becomes more effective since the mobility increases when the ECM decreases, and haptotaxis which tends to inhibit growth of the tumour into the less dense ECM outside the tumour (recall that haptotaxis induces motion up ECM gradients). Further, the MDE also degrades the pre-existing vessels which results in a reduction in the supply of oxygen. As a result of haptotaxis and the reduced oxygen supply, the tumour actually shrinks slightly after reaching a maximum radius of about 0.64mm, see Fig. 7(b). 4.2 Tumour-induced angiogenesis and vascular growth: No solid pressure-induced neo-vascular response We next consider tumour-induced angiogenesis where there is no effect of the solid pressure on either the radius of the neo-vessels or the extravasation of nutrient. In particular, we take  EMBED Equation.DSMT4 in Eqs. (4) and (27). Angiogenesis is initiated from the avascular tumour configuration at t=45days from Fig. 7. At this time, 10 sprout tips are released from the parent vessel. The initial vessel radii are 6mm. The inlet pressure and outlet pressures in the parent vessel are  EMBED Equation.DSMT4 and  EMBED Equation.DSMT4 respectively. The growth duration is t=0.05 which means that the intravascular flow and vessel adaption algorithms are called nearly every tumour growth time step. The flow duration is t=0.25 with a time step approximately equal to Dt=0.005 (again Dt is adaptive to satisfy an intravascular CFL condition). This means that 50 iterations of the flow and vascular adaptation algorithms are performed every tumour growth time step. By flowing and adapting the vascular network so frequently, we hoped that a relatively short flow duration time could be used to get a reasonable approximation of the blood flow in the network. Indeed, preliminary simulations showed that increasing the flow duration did not change the results qualitatively or, in some cases depending on the vascular network configuration, quantitatively. In a future work, we will quantify the effect of the flow duration upon the results. The evolution of the tumour and the neo-vascular network is shown in Figs. 8(a) and (b). As can be seen from the Figures, it takes some time for flow to develop after angiogenesis is initiated. Further, flow first occurs after about 7 days in a region near the parent vessel. This can be seen from the plots of haematocrit and oxygen which are signatures of blood flow. Little additional oxygen diffuses to the tumour. Accordingly, the tumour maintains a steady size (or shrinks a little due to the reasons described above). This may be seen in Fig. 8(c). Some of the vessels continue to lengthen, branch and migrate towards the tumour heading in particular for the hypoxic region where TAF is released. After about 10 days, a large loop forms through which blood flows. The loop penetrates the tumour and provides the tumour cells with a direct source of oxygen. The tumour responds by rapidly growing along the oxygen source and co-opts the neo-vessels and the hypoxic region shrinks and changes shape. As the tumour grows, the hypoxic and necrotic regions start to grow again as well and the neo-vessels near the tumour/host interface branch in response to wall shear stresses and increased TAF levels. This results in increased anastomosis and blood flow. The increased oxygen supply in turn causes large pressures to form in the proliferating region and the tumour to grow even more rapidly, enhancing this effect. Because there is no response of the neo-vessels to these large pressures, the tumour simply continues to co-opt the vessels creating an effective tumour microvasculature. This microvasculature provides a nearly uniform source of nutrient in the upper two thirds of the tumour; the lower third is primarily hypoxic and quiescent. As a consequence, the tumour shape remains compact as the tumour grows. In Fig. 8(b), the dimensional neo-vessel radii (in m) and intravascular pressures (in Pa) are shown together with the non-dimensional ECN and TAF concentrations. At early times, the radii are small and TAF diffuses from the quiescent zone. The ring of lowered ECM surrounding the tumour is clearly seen. The pressure is highest in the neo-vessels closest to the inlet of the parent capillary where the highest pressures are. As blood flow starts, the radii increase and the overall pressure decreases while the pressure in some vessels increases as blood spreads throughout the network. This process continues as the tumour grows and the vasculature continues to branch, anastomose and carry more and more flow. As the hypoxic and necrotic regions shrink, the TAF distribution changes and the vessels respond accordingly. Observe that the degraded ECM just outside the tumour does not prevent the vessels from penetrating the tumour even though the sprout-tips have to migrate up ECM gradients to accomplish this. The first vessels that penetrate the tumour do not carry blood and thus the tumour does not respond to their penetration. Instead, these vessels migrate towards the hypoxic region where they tend to get stuck. This occurs because the TAF concentration is nearly uniform (T=1) and so the sprout-tips to move randomly and tend to collide with their own trailing vessel preventing further migration. At later times though, new vessels grow into the tumour center and anastomose. This leads to blood flow and oxygen extravasation deep in the tumour interior. Further, observe that the tumour grows so fast that it outruns the ring of degraded ECM around its boundary and is growing into only very slightly degraded ECM. The ECM in the tumour interior degrades rather slowly and the ECM signature of the original avascular tumour spheroid can still be seen at late times. This simulation shows that when the neo-vessels are not affected by the tumour solid pressure, dramatic growth occurs as the tumour co-opts the host vasculature to create its own microvasculature and receives a direct source of oxygen. In addition, the tumour growth and angiogenesis processes are nonlinearly coupled as the vasculature responds to the growth by migrating towards the ever changing TAF distributions and by branching and anastomosing near the tumour-host interface. This leads to increased blood flow. At the same time, the increased blood flow in the vascular network affects how the tumour grows, and in particular speeds growth up. This then affects the response of the vasculature. 4.3 Tumour-induced angiogenesis and vascular growth: The effect of solid pressure-induced neo-vascular response Next, we consider, in Figs. 9(a)-(c), the effect of solid/mechanical pressure-induced vascular response on tumour-induced angiogenesis and vascular growth. We repeat the simulation in section 4.2 except with  EMBED Equation.DSMT4 non-zero as given in Eq. (5). This means that transfer of oxygen from the neo-vessels to the tissue may be significantly reduced and the vessel radii may be correspondingly constricted. With the values of the parameters used here (Tables 4-6), a solid pressure-induced vascular response begins to occur when the solid pressure EMBED Equation.DSMT4 . At early times, the angiogenic response and the tumour growth is similar to the case presented earlier in Figs. 8(a)-(c). The newly developing vessels migrate, proliferate, branch and anastomose. It also takes some time for flow to begin with significant flow developing only after about 10 days. Blood flow in the neo-vasculature starts near the parent capillary and eventually the flow reaches the tumour. Because the initial ECM is slightly different than that in Fig. 8 (due to the random component) and due to the random component of the sprout tip motion, the vascular network at early times is not identical to that obtained previously in Fig. 8. In contrast to the case considered in Fig. 8, here the solid pressure prevents any delivery of oxygen internally to the tumour and thus the delivery of oxygen is heterogeneous and significant oxygen gradients persist in the tumour interior. There is no functional microvasculature internal to the tumour. While the tumour responds by growing towards the oxygen-delivering neo-vessels, the solid pressure generated by tumour cell proliferation also constricts the neo-vessels in the direction of growth (where pressure is highest) and also correspondingly inhibits the transfer of oxygen from those vessels. As a consequence, the overall solid pressure is significantly lower than that in Fig. 8. This makes the tumour grow much more slowly than that in Fig. 8 as can be seen in Fig. 9(c). Note the vertical scale in Fig. 9(c) is one half of that in Fig. 8(c). The neo-vessels in other areas of the host microenvironment then provide a stronger source of oxygen. This triggers tumour-cell proliferation and growth in regions where proliferation had been decreased previously. The heterogeneity of oxygen delivery and the associated oxygen gradients cause heterogeneous tumour cell proliferation. Unlike the case in Fig. 8, proliferation is confined to regions close to the tumour-host interface. This results in morphological instability that leads to the formation of invasive tumour clusters (e.g. buds) and a complex tumour morphology. This result is consistent with the theory and predictions made earlier (see, for example, Cristini et al. (2003), Cristini et al. (2005), Anderson (2005), Macklin and Lowengrub (2005, 2006, 2007), Anderson et. al. (2006) and Gerlee and Anderson (2007)), that substrate inhomogeneities in the tumour microenvironment tend to cause morphological instabilities in growing tumours. Although nutrient-providing, functional vessels are not able to penetrate the tumour during growth, the growth of the tumour elicits a strong branching and anastomosis response from nearby neo-vessels in the host microenvironment. Although there is an analogous neo-vascular response seen in Fig. 8, the effect here is much more pronounced as the levels of TAF are higher in these regions (because tumour hypoxia is increased) and thus the wall shear stresses initiate more significant branching. In Fig. 9(b), the dimensional neo-vessel radii (in m) and intravascular pressures (in Pa) are shown together with the non-dimensional ECN and TAF concentrations. As before, blood flow causes a dilation of the vessels and an overall decrease of pressure as branching, anastomosis and increased blood flow occurs throughout the neo-vascular network. The constriction of neo-vessels in response to the solid pressure is clearly seen. The tumour-secreted MDE degrades the ECM in the host microenvironment near the tumour and in the tumour interior. As before (recall Fig. 8(b)), the neo-vessels are still able to migrate through the region of lower ECM even though this acts against haptotaxis. Because the tumour grows more slowly than that in Fig. 8(b), only the tips of the invasive clusters outrun the degraded ECM. As can be seen in Fig. 9(c), the host ECM is degraded in the region between the invading clusters. The ECM signature of the original avascular tumour spheroid can no longer be seen at later times. This simulation shows even stronger nonlinear coupling between the tumour-induced angiogenesis and the progression of the tumour compared to the prior case shown in Fig. 8. The pressure-induced vascular response of constricting neo-vessel radii and inhibiting blood-tissue oxygen transfer not only affects the tumour growth dramatically, but also significantly affects the growth of the neo-vascular network, and vice-versa. 4.4 Avascular growth to a multicellular (2D) spheroid with enhanced ECM degradation and production We next examine the effect of ECM degradation upon the results. In Fig. 10(a), we repeat the simulation in section 4.1 except that both the MDE degradation and production parameters are increased (see Table 4). The tumour grows by uptaking oxygen delivered by the pre-existing (uniform) vasculature and growth is more rapid than that for the avascular tumour shown in section 4.1 (Fig. 7(a)). This occurs because the mobility is larger here due to the enhanced degradation of ECM. This effect overcomes the tendency of haptotaxis to keep the tumour away from the degraded ECM. The tumour reaches a nearly steady size, containing both a hypoxic and a necrotic core, that is significantly larger than that shown in Figs. 7(a) and (b); the radius at 45 days is approximately 0.78 mm (see Fig. 10(b)). At the final time shown (45 days), the ECM is significantly degraded in the host microenvironment and in the tumour necrotic core to the point that there is even a thin annular hole in the ECM, immediately surrounding the spheroid, and a circular hole in the necrotic region, where the density of ECM  EMBED Equation.DSMT4 . 4.5 Tumour-induced angiogenesis and vascular growth: The effect of solid pressure-induced neo-vascular response and enhanced ECM degradation and production Next, we consider, in Figs. 11(a)-(c), the effect of enhanced ECM degradation on tumour-induced angiogenesis and vascular growth. We repeat the simulation shown in section 4.3 except that the initial condition is the t=45 day simulation from Fig. 10(a) and the MDE parameters are the same as in section 4.4 (see Tables 4-6). As in the simulation shown in section 4.1, the neo-vessels grow and form loops near the parent capillary. However, now because of the growing ECM annular hole surrounding the tumour, the neo-vessels are not able to reach the tumour and are instead trapped by the ECM hole due to haptotaxis. The vessels then encapsulate roughly the upper half of the tumour. As blood flows through the neo-vascular network and approaches the tumour, the tumour responds by growing towards the flowing neo-vessels that provide the oxygen source, as in section 4.3. The tumour elongates, constricts the neo-vessels in its path and prevents the transfer of oxygen from the neo-vessels to the host. This limits tumour cell proliferation and results in a roughly steady maximum solid pressure. Correspondingly, there is heterogeneous oxygen supply, heterogeneous tumour cell proliferation and there are strong oxygen gradients. As in section 4.3, this results in a morphological instability of the growing solid tumour. As the tumour continues to grow, the neo-vessels respond by increasing branching and anastomosing near the tumour-host interface. This results in a broader supply of oxygen in the part of the tumour closest to the parent capillary. Proliferation is increased and the top of the tumour flattens. The increased proliferation leads to large solid pressures which then constrict the nearby neo-vessels and inhibit oxygen supply. The tumour then begins to grow towards other vessels near the parent capillary and the top of the tumour becomes unstable. Further, there is instability along the side of the tumour that leads to the encapsulation of host domain inside the tumour. Also observe that a small amount of oxygen is able to be delivered to the tumour at very late times as haematocrit is trapped in a constricted vessel at a location where the pressure is sufficiently low to allow extravasation. Fig. 11(b) shows the dimensional neo-vessel radii (in m) and the intravascular pressures (in Pa) together with the non-dimensional ECN and TAF concentrations. The results are similar to those obtained before except that the tumour does not outrun the ECM hole although at the top of the tumour, the hole is quite shallow. Interestingly, even though the initial tumour in Fig. 11(a) is larger than that in Fig. 9(a), the final tumour size at t=150 days is roughly the same for both cases (see Figs. 11(c) and 9(c)). The ECM hole present in the simulation in Fig. 11 prevents the neo-vessels from getting close to the tumour during the early stages of growth; this allows the tumour in Fig. 9 to catch up and even grow slightly larger than that in Fig. 11. Finally, in Fig. 12, we compare the average radii in the neo-vascular networks for the simulations in Figs. 8, 9 and 11. At early times, the radii for the simulation in Fig. 8, where the neo-vasculature does not respond to solid pressure, grows the fastest as blood flows uninhibited through the network. Later, however, the simulation with lower ECM degradation shows the most rapid radii increase. This occurs because the EC sprout-tips are able to move more freely through the host domain and do not get caught by degraded ECM. This provides the vascular network with a more widely varying flow response. 5. Conclusions and future directions In this paper we have coupled an improved continuum model of solid tumour invasion (following Macklin and Lowengrub 2007) with a model of tumour-induced angiogenesis (following McDougall et al. (2006)) to produce a new multi-scale model of vascular solid tumour growth. The invasion and angiogenesis models were coupled through the tumor angiogenic factors (TAF) released by the tumor cells and through the nutrient extravasated from the neo-vascular network. As the blood flows through the neo-vascular network, nutrients (e.g. oxygen) are extravasated and diffuse through the ECM triggering further growth of the tumour, which in turn influences the TAF expression. In addition, the extravasation is mediated by the hydrostatic stress (solid pressure) generated by the growing tumour. The solid pressure also affects vascular remodelling by restricting the radii of the vessels and thus the flow pattern and wall shear stresses. The vascular network and tumour progression were also coupled via the ECM as both the tumour cells and ECs upregulate matrix degrading proteolytic enzymes which cause localized degradation of the ECM which in turn affects haptotactic migration. We performed simulations of the multi-scale model that demonstrated the importance of the nonlinear coupling between the growth and remodelling of the vascular network, the blood flow through the network and the tumour progression. The solid pressure generated by tumour cell proliferation effectively shuts down large portions of the vascular network dramatically affecting the flow, the subsequent network remodelling, the delivery of nutrients to the tumour and the subsequent tumour progression. In addition, ECM degradation by tumour cells was seen to have a dramatic affect on both the development of the vascular network and the growth response of the tumour. In particular, when the ECM degradation is significant, the newly formed vessels tended to encapsulate, rather than penetrate, the tumour and were thus less effective in delivering nutrients. There are many directions in which this work will be taken in the future both in terms of modelling additional biophysical effects as well as algorithmic improvements. Regarding the algorithm, we plan to upgrade the solid pressure/nutrient solver by solving for P and s as a coupled system. This will prevent oscillations that may occur by lagging P in the source term for nutrient. We also plan to accelerate the solver for the intravascular pressure to improve performance of the coupled algorithm. Regarding the model, we plan to develop a more detailed model of the effect of solid pressure on the constriction and collapse of vessels in the microvasculature and on the corresponding response of the microvascular network. We also plan to include, for the first time, the effects of the venous system. Other features, such as the recruitment of pericytes by the vascular ECs will also be investigated. In addition, we will incorporate more realistic models for soft tissue mechanics. The work presented here demonstrates that nonlinear simulations are a powerful tool for understanding phenomena fundamental to solid tumour growth. A biophysically justified computer model could provide an enormous benefit to the clinician, the patient, and society by efficiently searching parameter space to identify optimal, or nearly optimal, individualized treatment strategies involving, for example, chemotherapy and adjuvant treatments such as anti-angiogenic or anti-invasive therapies. This is a direction we plan to explore in the future. Acknowledgements ML acknowledges partial support from a U.S. Department of Education GAANN (Graduate Assistance in Areas of National Need) fellowship. JSL acknowledges partial support from the National Science Foundation (NSF) through grants DMS-0352143 and DMS-0612878 and from the National Institutes of Health (NIH) grant P50GM76516 for a National Center of Excellence in Systems Biology at UCI. VC acknowledges partial support from the NSF and the (NIH) through grants NSF-DMS-0314463 and NIH-5R01CA093650-03. Further, ML, VC and JL are grateful to Herman Frieboes for valuable discussions. Appendix 1: Details of the stimuli that affect the vessel radius in the adaptive dynamic tumour-induced angiogenesis model Wall shear stress Many studies show that vessels adapt their radius in order to maintain a constant level of wall shear stress (e.g., Pries et al. (1998, 2001a, 2001b), Fung 1993) . Hence vessel radius tends to increase with increasing wall shear stress, whilst wall shear stress decreases with increasing radius. The non-dimensional wall shear stress stimulus can be described by a logarithmic law as  EMBED Equation.DSMT4 , (A1.1) where  EMBED Equation.DSMT4  is the actual wall shear stress in a vessel segment,  EMBED Equation.DSMT4 is the flow rate in the vessel under consideration, tref is a constant included to avoid singular behaviour at low shear rates, and EMBED Equation.DSMT4  is a wall shear stress scale factor (Pries 2001a). Dimensional stresses are in dynes/cm2.The dimensional wall shear stress calculated in the parent vessel of our computational model is of the order 4Pa (40dynes/cm2) and capillary values are less than 2Pa (20dynes/cm2), in agreement with those measured experimentally in the dog by Kamiya et al. (1984). Adaptation in response to the wall shear stress stimulus alone tends to reinforce a single path in the network composed of a few well-established fully dilated vessels corresponding to the main flowing backbone of the vasculature whilst simultaneously eliminating the low-flow paths. However, the resulting network is unstable in the sense that there is no consistent balance for the radius and flow distribution achieved when Swss is considered in isolation. Intravascular pressure Intravascular pressure is another key stimulus for vascular adaptation. Pries et al. (1995) have experimentally observed on the rat mesentery the dependence of the magnitude of the wall shear stress with the local intravascular pressure (Pvessel). They proposed a parametric description of their experimental data, which exhibits a sigmoidal increase of the wall shear stress with increasing pressure through the following. The sensitivity of the corresponding (non-dimensionall) stimulus to intravascular pressure is then described by:  EMBED Equation.DSMT4 , (A1.2) where  EMBED Equation.DSMT4 , EMBED Equation.DSMT4 is a intravascular pressure stress scale factor (note that  EMBED Equation.DSMT4  itself is not needed for the simulation as the above formula is used),  EMBED Equation.DSMT4 is a pressure scale factor, and kp is rate of response of the radius to the pressure stimulus. Metabolic haematocrit-related stimulus The metabolic stimulus effectively stabilises the adapting network by stimulating vessel growth in areas of the vascular bed exhibiting low flow. The non-dimensional stimulus is once again described by a logarithmic law given by:  EMBED Equation.DSMT4 , (A1.3) where km is a constant characterizing the relative intensity of the metabolic stimulus, and  EMBED Equation.DSMT4 is the flow rate in the parent vessel. Natural shrinking tendency The natural reaction of the basal lamina is thought to counter any increase in vessel diameter and can be modelled by (Pries et al. 1998):  EMBED Equation.DSMT4 , (A1.4) where EMBED Equation.DSMT4 ,  EMBED Equation.DSMT4 is a constant and  EMBED Equation.DSMT4 is a shrinking tendency scale factor (taken to be  EMBED Equation.DSMT4 Pa). Note that since  EMBED Equation.DSMT4 =1 Pa,  EMBED Equation.DSMT4 no longer appears in Eq. (26) and instead we obtain EMBED Equation.DSMT4 . To summarise, the conditions for vessel branching in the DATIA model are as follows: (i) the likelihood of a vessel branching increases with both the local TAF concentration and the magnitude of the shear stress affecting the vessel wall; (ii) the vessel must reach a certain level of maturation before it is able to branch, although branching cannot occur once a basal lamina has formed around a vessel (Gee et al. 2003, Benjamin et al. 1998). In the simulations presented here, however, there is no upper limit on the age of vessels that are allowed to branch. This models the case in which the continued tumour-induced angiogenic response prevents the formation of basal lamina. This effect will be considered in a future work. Appendix 2: Numerical techniques for tumour invasion model The ghost cell method The reaction-diffusion equations describing nutrient transport, the pressure equilibration, the matrix degrading enzyme evolution and tumour angiogenic factor distribution may be written generically in the form:  EMBED Equation.DSMT4 , (A2.1) where v = s, P, M or T and t is either a pseudo-time variable (in the case of the elliptic/quasi-steady equations for nutrient, pressure and angiogenic factor) or is the real time (in the case of the MDE equation). Note that without loss of generality, fR is assumed to be non-positive. We discretize Eq. (A2.1) in (pseudo/real) time using the backward Euler method and lag the dependence of the diffusion coefficient D, reaction fR and source fS terms upon v. Away from interfaces, Eq. (A2.1) is discretized in space using a centered finite difference scheme. For quantities that are smooth (i.e. s, M, E or T) centered differences are used throughout the domain. For the pressure, however, the jump conditions in Eqs. (11) and (12) need to be correctly incorporated and the difference stencils are modified near the  accordingly (Macklin and Lowengrub 2006, 2007). Note that we assume that there are no jumps across N and thus no stencil modification is needed near N. The stencil is modified in the following way. Whenever the boundary  intersects the computational stencil, we replace points that are outside of the boundary with  ghost points (denoted by hats) that are extrapolated from within the region where we also use the jump boundary condition, i.e.  EMBED Equation.DSMT4 , (A2.2) where the  sgn() ensures the jump condition is applied in the proper direction and is discretized using the approach in (Sussman & Fatemi 1999). See Figure 3. In Eq. (A2.2), the curvature k is accurately obtained from the level-set function  (i.e.,  EMBED Equation.DSMT4 ) even for complex interface morphologies using a geometry aware discretization we recently developed (Macklin & Lowengrub 2006, 2007). Further, Eq. (A2.2) introduces non-grid points Pr and Pl into the discretization. However, taking into account the jump condition on the normal derivative in Eq. (12), one obtains two equations for the two variables Pl and Pr allowing them to be eliminated from the discretization. The normal vector  EMBED Equation.DSMT4 is also discretized using a geometry-aware method (Macklin & Lowengrub 2007). The normal jump in Eq. (12) is discretized by taking decomposing the normal derivative into components in the grid (ul, ur) and off-grid (diagonal, vl, vr) directions (see Figure 4) that involve pl and pr and carefully match extrapolations from inside and outside of the tumor domain . This approach does not smear tangential derivatives a problem that plagued earlier implementations of ghost cell methods. Indeed, this algorithm has tested better than 1.5-order accurate (Macklin & Lowengrub 2007); in contrast, previous ghost fluid and ghost cell methods (e.g., Liu et al. 2000) smear jumps in the tangential derivative and either fail to converge or are at best very low (sub 0.3) order accurate. The fully discrete version of Eq. (A2.1) is solved using a nonlinear adaptive Gauss-Seidel-type iterative method (NAGSI). Because of the lagged form of the time discretization, there is no need to solve a large linear system at each iteration step. Instead, there is only a local solve which can be performed analytically. In the case of the time-dependent MDE equation, this is sufficient to advance to the next time step. For the nutrient, pressure and tumour angiogenic factor, however, Eq. (A2.1) must be solved to steady state (within an error tolerance) and thus further iterations are necessary. To advance to steady state efficiently, we take advantage of the fact that as the solution converges, the numerical solution tends to change most on a small fraction of the computational nodes. Therefore, we can select computational nodes where the numerical solution is changing most rapidly and update only those nodes. This algorithm is very efficient and adaptivity typically yields about a 50% reduction in computational time (Macklin & Lowengrub 2007). The level-set method Level set methods were first developed by Osher and Sethian (Osher & Sethian 1988) and have been used to study the evolution of moving, complex surfaces (see the books (Sethian 1999, Osher & Fedkiw 2002) and the references therein). Level set methods have been previously used in the context of tumor growth by Macklin & Lowengrub (2005, 2006, 2007), Zheng et al. (2005), and Hogea et al. (2006). In the level set approach, instead of explicitly tracking the position of interfaces  and N and manually handling topology changes, the level set function is updated by solving a Hamilton-Jacobi partial differential equation, which automatically accounts for the interface motion and possible morphological changes of the interface. The idea is as follows. Let V be the outward normal velocity of the interface S (e.g. see Eq. (13)) and let Vextension be an extension of V off of the tumor boundary  such that Vextension = V on . One possibility is to use u, from Eq. (7), as the extension of V off the interface. This is actually not done, as explained below. To update the position of the boundary , we solve the Hamilton-Jacobi equation  EMBED Equation.3 . (A2.3) Note that the position of the necrotic/viable interface N is similarly obtained by evolving N with an extension VN,extension of the necrotic velocity VN off the boundary N. We obtain the extension velocities using a bilinear extrapolation technique first presented in Macklin and Lowengrub (2005) that is constant in the normal direction to the interface. This can be thought of as a higher order version of the fast marching method (e.g. Adalsteinsson and Sethian 1999, Zhao 2005). An advantage of taking Vextension to be constant in the normal direction is that distance functions are then preserved by Eq. (A2.3). If instead Vextension = u, this would not be the case. Even though theoretically  EMBED Equation.DSMT4 , it is still helpful to prevent accumulation of numerical error by re-initializing  periodically to be a distance function by solving  EMBED Equation.DSMT4 , (A2.4) to steady state, where t is pseudo-time and  EMBED Equation.DSMT4 is the original level-set function prior to reinitialization. Finally, Eqs. (A2.3) and (A2.4) are discretized in time using first order Euler method and a third-order total variation-diminishing Runge-Kutta method (TVD-RK; Gottleib & Shu 1997, Gottlieb et al. 2001), respectively. The former is used to reduce computational cost. In space, the  EMBED Equation.DSMT4  is discretized using the fifth-order weighted essentially non-oscillatory (WENO) method (Jiang & Shu 1996, Jiang & Peng 2000). Lastly, we update Eqs. (A2.3) and (A2.4) only near the interface using the narrow band/local level set technique (Malladi et al. 1996, Peng et al. 1999) with a band size of 20Dx. For further details on the ghost cell/level-set algorithm, the reader is referred to Macklin and Lowengrub 2005, 2006 and 2007. Appendix 3: Additional details of the invasion model: The interpolating functions Several interpolation functions are used in the tumour invasion model. These are given below. In the oxygen uptake equation (2), the interpolation function ps(E0) is defined by  EMBED Equation.DSMT4 , and p(E0) is a cubic polynomial satisfying:  EMBED Equation.3 . In Eq. (2), q() is a quartic polynomial chosen to satisfy:  EMBED Equation.3  In Eq. (5), pcutoff is a cubic polynomial chosen to satisfy:  EMBED Equation.3 . Both the haptotactic sensitivity and the mobility of the tumour cells involve interpolating functions. In Eq. (15), pc (E) is a quartic polynomial chosen to satisfy:  EMBED Equation.3  Finally, in Eq. (16), the function pm (E) is the cubic polynomial:  EMBED Equation.3  References V. Cristini, J. Lowengrub and Q. Nie. Nonlinear simulation of tumour growth. J. Math. Biol. 46(2003), 191-224. Zheng, X., Wise, S.M., Cristini, V. Nonlinear simulation of tumor necrosis, neo-vascularization and tissue invasion via an adaptive finite-element/level-set method (2005) Bull. Math. Biol. 67: 211-259 Araujo R.P. and McElwain D.L.S. 2004. A history of the study of solid tumor growth: the contribution of mathematical modelling. Bull. Math. Biol. 66:1039-1091. Araujo R.P., McElwain D.L.S. 2005. A mixture theory for the genesis of residual stresses in growing tissues I: A general formulation. SIAM J. Appl. Math. 65:1261-1284. Araujo R.P., McElwain D.L.S. 2005. A mixture theory for the genesis of residual stresses in growing tissues II: Solutions to the biphasic equations for a multicell spheroid. SIAM J. Appl. Math. 66:447-467. Deakin AS. 1976. Model for initial vascular patterns in melanoma transplants. Growth 40: 191-201. Hanahan D., Weinberg R.A. 2000. The hallmarks of cancer. Cell 100:57-70. Ambrosi D., Preziosi L. 2002. On the closure of mass balance models for tumor growth. Math Mod. Meth. Appl. Sci. 12:737-754. Byrne H.M., Preziosi L. 2003. Modeling of solid tumour growth using the theory of mixtures. Math. Med. Biol. 20:341:366. Balding D, McElwain DLS. 1985. A mathematical model of tumour-induced capillary growth. J. theor. Biol. 114: 53-73. Stokes CL, Lauffenburger DA. 1991. Analysis of the roles of microvessel endothelial cell random motility and chemotaxis in angiogenesis. J. Theor. Biol. 152: 377-403. Friedl P., Wolf K. 2003. Tumor-cell invasion and migration: Diversity and escape mechanisms. Nature Reviews 3:362-371. Carmeliet P. 2005. Angiogenesis in life, disease and medicine. Nature 438:932-936. Chaplain MAJ, Stuart AM. 1993. A model mechanism for the chemotactic response of endothelial cells to tumour angiogenesis factor. IMA J. Math. Appl. Med. Biol. 10: 149-68. Fung Y.C. 1993 Biomechanics. New York, Springer. Chomyak O.G., Sidorenko M.V. 2001. Multicellular spheroids model in oncology. Exp. Oncology 23:236-241. Kunz-Schughard L.A., Freyer J.P., Hofstaedter F., Ebner R. 2004. The use of 3-D cultures for high-throughput screening: The multicellular spheroid model. J. Biomol. Screening 9:273-285. Kim J.B. 2005. Three-dimensional tissue culture models in cancer biology. Seminars in Cancer Biology 15:365-377. Walles T., Weimer M., Linke K., Michaelis J., Mertsching H. 2007. The potential of bioartificial tissues in oncology research and treatment. Onkologie 30:388-394. Byrne HM, Chaplain MAJ. 1995. Mathematical models for tumour angiogenesis: numerical simulations and nonlinear wave solutions. Bull. Math. Biol. 57: 461-86. Orme ME, Chaplain MAJ. 1997. Two-dimensional models of tumour angiogenesis and anti-angiogenesis strategies. IMA J. Math. Appl. Med. Biol. 14: 189205. Adalsteinsson D., Sethian J.A. 1999. The fast construction of extension velocities in level-set methods. J. Comput. Phys. 148:2-22. Olsen L, Sherratt JA, Maini PK, Arnold F. 1997. A mathematical model for the capillary endothelial cell-extracellular matrix interactions in wound-healing angiogenesis. IMA J. Math. Appl. Med. Biol. 14: 26181. Anderson ARA, Chaplain MAJ. 1998. Continuous and discrete mathematical models of tumor-induced angiogenesis, Bull. Math. Biol. 60: 857-99. Stephanou A., McDougall S.R., Anderson A.R.A. and Chaplain M.A.J. 2006. Math. Comp. Model. 44:96-123. Palecek SP, Loftus JC, Ginsberg MH, Lauffenburger DA, Horwitz AF. 1997. Integrin-ligand binding properties govern cell migration speed through cell-substratum adhesiveness. Nature 385:537-540. DiMilla P.A., Barbee K., Lauffenburger D.A. 1991 Mathematical model for the effects of adhesion and mechanics on cell migration speed, Biophys. J. 60:15-37. Dickinson R.B. and Tranquillo R.T. 1993 A stochastic model for adhesion-mediated cell random motility and haptotaxis, J. Math. Biol. 31:563-600. Chaplain MAJ. 2000. Mathematical modelling of angiogenesis. Journal of Neuro-Oncology 50: 37-51. Sanga S., Sinek J., Frieboes H., Ferrari M., Fruehauf J., Cristini V. 2006. Mathematical modelling of cancer progression and response to chemotherapy Levine HA, Pamuk S, Sleeman BD, Nielsen-Hamilton M. 2001. Mathematical modeling of the capillary formation and development in tumor angiogenesis: penetration into the stroma. Bull. Math. Biol. 63: 801-63. Kloner R.A., Jennings R.B. 2001. Consequences of brief ischemia: stunning, preconditioning and their clinical implicatlins: Part 1. Circulation 104:2981-2989. Galaris D., Barbouti A., Korantzopoulos P. 2006. Oxidative stress in hepatic ischemiareperfusion injury: the role of antioxidants and iron chelating compounts. Current Pharma. Design 12:2875-2890. Plank MJ, Sleeman BD. 2004. Lattice and non-lattice models of tumour angiogenesis. Bull. Math. Biol. 66: 1785-1819. Zhao H.K. 2005. A fast sweeping method for Eikonal equations. Math. Comp. 74:603-627. Byrne H.M., Alarcon T., Owen M.R., Webb S.D. and Maini P.K. 2006. Modeling aspects of cancer dynamics: A review. Phil. Trans. R. Soc. A 364:1563-1578. Graziano L., and Preziosi L. 2007. Mechanics in tumor growth, in Modelling of Biological Materials, ed. F. Mollica, K.R. Rajagopal and L. Preziosi, Birkhauser, 267-328. Alarcon T., Byrne H.M., Maini P.K.. 2005 A multiple scale model for tumor growth, Multiscale Mod. Simul. 3:440-475. Quaranta V., Weaver A.M., Cummings P.T., and Anderson A.R.A. 2005. Mathematical modelling of cancer: The future of prognosis and treatment, Clin. Chim. Acta 357:173-179. Gerlee P., Anderson A.R.A. 2007. Stability analysis of a hybrid cellular automaton model of cell colony growth. Phys. Rev. E 75:0151911. Roose T., Chapman S.J., and Maini P.K. 2007. Mathematical models of avascular cancer. SIAM Reviews. 49:179-208. Mantzaris NV, Webb S, Othmer HG. 2004. Mathematical modeling of tumor-induced angiogenesis. J. Math. Biol. 49: 111-87. L.A. Liotta, G.M. Saidel and J. Kleinerman, Diffusion model of tumor vascularization and growth. Bull. Math. Biol. 39: 117-129 (1977) H.P.Greenspan, On the growth and stability of cell cultures and solid tumours. J. theor. Biol. 56: 229-242 (1976). M.E. Orme \& M.A.J. Chaplain A Mathematical Model of Vascular Tumour Growth and Invasion (1996) Math. Comp. Modelling 23, 43-60. T. Alarcon, H.M. Byrne, P.K. Maini. A cellular automaton model for tumour growth in a heterogeneous environment. J. theor. Biol. (2003) 225, 257274. Kaur B., Khwaja F.W., Severson E.A., Matheny S.L., Brat D.J., VanMeir E.G. 2005. Hypoxia and the hypoxia-inducible-factor pathway in glioma growth and angiogenesis. Neuro-Oncol. 7:134-153. Erler J.T., Bennewith K.L., Nicolau M., Dornhoefer N., Kong C.L., Chi Q-T., Jeffrey S.S., Siaccia A.J. 2006. Lysyl oxidase is essential for hypoxia-induced metastasis. Nature 440:1222-1226. Pouyssegur J., Dayan F., Mazure N.M. 2006. Hypoxia signalling in cancer and approaches to enforce tumour regression. Nature 441:437-443. Some more references: Alarcon, T., Byrne, H., Maini, P., 2003. A cellular automaton model for tumour growth in inhomogeneous environment. J. Theor. Biol. 225 (2), 257274. Anderson, A.R.A., 2005. A hybrid mathematical model of solid tumour invasion: the importance of cell adhesion. Math. Med. Biol. 22, 163186. Anderson, A.R.A., Chaplain, M.A.J., Newman, E.L., Steele, R.J.C., Thompson, A.M., 2000. Mathematical modelling of tumour invasion and metastasis. J. Theor. Med. 2, 129154. Anderson, A.R.A., Chaplain, M.A.J., 1998. Continuous and discrete mathematical models of tumor-induced angiogenesis. Bull. Math. Biol. 60, 857899. Ausprunk, D.H., Folkman, J., 1977. Migration and proliferation of endothelial cells in preformed and newly formed blood vessels during tumour angiogenesis. Microvasc. Res. 14, 5365. Folkman, J., 1971. Tumor angiogenesis: therapeutic implications. N. Eng. J. Med. 285, 11821186. Folkman, J., 1995. Angiogenesis in cancer, vascular, rheumatoid and other disease. Nature Med. 1, 2131. Folkman, J., Klagsbrun, M., 1987. Angiogenic factors. Science 235, 442447. Gimbrone, M.A., Cotran, R.S., Leapman, S.B., Folkman, J., 1974. Tumor growth and neovascularization: an experimental model using the rabbit cornea. J. Natl. Cancer Inst. 52, 413427. Godde, R., Kurz, H., 2001. Structural and biophysical simulation of angiogenesis and vascular remodeling. Dev. Dynam. 220, 387401. Jain, R.K., 1987. Transport of molecules in the tumor interstitium: a review. Cancer Res. 47 (12), 30393051. Jain, R.K., 1988. Determinants of tumour blood flow: a review. Cancer Res. 48 (10), 26412658. Jain, R.K., 2001. Normalizing tumor vasculature with anti-angiogenic therapy: a new paradigm for combination therapy. Nature Med. 7, 987989. Jain, R.K., 2003. Molecular regulation of vessel maturation. Nature Med. 9, 685693. Jain, R.K., 2004. Recognition of tumor blood vessel normalization as a new antiangiogenic concept. Nature Med. 10, 329330. Kamiya, A., Bukhari, R., Togawa, T., 1984. Adaptive regulation of wall shear stress optimizing vascular tree function. Bull. Math. Biol. 46, 127137. Krenz, G.S., Dawson, C.A., 2002. Vessel distensibility and flow distribution in vascular trees. J. Math. Biol. 44, 360374. Levine, H.A., Pamuk, S., Sleeman, B.D., Nielsen-Hamilton, M., 2001. Mathematical modeling of the capillary formation and development in tumor angiogenesis: penetration into the stroma. Bull. Math. Biol. 63, 801863. Mantzaris, N.V., Webb, S., Othmer, H.G., 2004. Mathematical modelling of tumor-induced angiogenesis. J. Math. Biol. 49, 111187. McDougall, S.R., Sorbie, K., 1997. The application of network modelling techniques to multiphase flow in porous media. Petrol. Geosci. 3, 161169. McDougall, S.R., Anderson, A.R.A., Chaplain, M.A.J., Sherratt, J.A., 2002. Mathematical modelling of flow through vascular networks: implications for tumour-induced angiogenesis and chemotherapy strategies. Bull. Math. Biol. 64 (4), 673702. Morikawa, S., Baluk, P., Kaidoh, T., Haskell, A., Jain, R.K., McDonald, D.M., 2002. Abnormalities in pericytes on blood vessels and endothelial sprouts in tumors. Am. J. Pathol. 160, 9851000. Muthukkaruppan, V.R., Kubai, L., Auerbach, R., 1982. Tumor-induced neovascularization in the mouse eye. J. Natl. Cancer Inst. 69, 699705. Olsen, L., Sherratt, J.A., Maini, P.K., Arnold, F., 1997. A mathematical model for the capillary endothelial cell-extracellular matrix interactions in wound-healing angiogenesis. IMA J. Math. Appl. Med. Biol. 14, 261281. Orme, M.E., Chaplain, M.A.J., 1997. Two-dimensional models of tumour angiogenesis and anti-angiogenesis strategies. IMA J. Math. Appl. Med. Biol. 14, 189205.  #.:GIJOVWYajkln}úð˥˜ˑÜ|tfh h 6H*mH sH hMmH sH h hMH*mH sH hH*mH sH hM^BhDmH sH h H*mH sH hM^BhMmH sH hv0JmH sH h#H*mH sH h6] mH sH h#mH sH hM^Bhw mH sH h hM^Bhh}EhdL hI hM^Bhw hl$IJ x D gxzb%c%gdMiE$a$gd0Qugd,_$a$gdS$a$gdw $ a$gdygd6] gdw gdw % '&''  x y 7 8 9 C D E  @  % 3 9 A N Q [ */Pż{vnianai{i{ih]hw 5 h45h]h&5 he15 hy5h]hMiE56 h 5h]hMiE5h hM^Bhw h6] 5\mH nHsH tHh hmH sH h6mH sH hh6H*mH sH h 6mH sH hU6mH sH h h 6H*mH sH h h 6mH sH %6IJV_m-m:^`dy )7f½̸̳|un hvhB hvhw hvh&hvhw 5hvh]5hM^Bhw 5aJhM^Bhw aJh]hy5 h{5 hy5 hr5 hzV5 hoc5 h6q5 h$5 h 5 hF5 h5 he15h]hw 5h]h&5 h45 h.5*fghkxyz{ @ef'7>hijxĵ||||tltd\h3kLmH sH hnFmH sH h8mH sH h,mH sH hmH sH hySmH sH hamH sH hM^Bh6mH sH hM^Bh0QumH sH hM^BhMiEmH sH hM^Bh0QuCJaJmH sH *hM^Bh0Qu5CJ\aJmH nHsH tH$hc'5CJ\aJmH nHsH tHh,_mH sH hShw CJaJ!7<=/67z{>?ZmqrIJWⶪⶒגhx)mH sH hPmH sH hM^BhUI6mH sH hM^Bh&6mH sH hM^BhMiE6mH sH hM^BhUImH sH hM^Bh&mH sH hM^Bh6mH sH hM^Bh0QumH sH hM^BhMiEmH sH hM^Bh8tmH sH hMiEmH sH 2W:;}%1]lz  BU⸩vvjh mH nHsH tHh*vmH nHsH tHhUohmH nHsH tHhb{QmH nHsH tHhM^Bh0QumH nHsH tHhM^Bh>mH nHsH tHhM^BhUImH nHsH tHhcmH sH hPmH sH hM^Bh0QumH sH hM^BhUImH sH hM^BhMiEmH sH hb{QmH sH &Pa # Z ` a b i q r !\!m!{!!!!!""# #!#*#.#5#˸˸˸˸}umeh1qmH sH hz!.mH sH hmH sH hY8mH sH hkmH nHsH tHhM^Bh?mH nHsH tHhVmH nHsH tHh mH nHsH tHh mH sH h h mH sH h hb{QmH sH hM^BhUImH nHsH tHhM^Bh0QumH nHsH tHhmH nHsH tH'5#?#Q#g#m#####"$#$h$i$$$4%5%a%b%c%%%%%a&j&&&&&'"'-'N''(7(W(X(Y(~(((() )c)d))))))Ǽh)mH sH hmH sH hDmH sH h0[hmH sH hmH sH hQhmH sH hM^BhmH sH hMiEmH sH hM^Bh0QumH sH hM^Bh>mH sH hM^BhMiEmH sH hmH sH hz!.mH sH 4)))))*8*d*+8+9+:+_+f+p+q+++1,2,T,,,,^-~-...W.X...../"/$/*/+/;/[/]/////I1J1K1L1·ʷʷʷʷʷʷhM^BhQmH sH homH sH h\mH sH h_mH sH hQh^tmH sH hQmH sH hQhQmH sH h0[mH sH h)mH sH hDmH sH h0[h0[mH sH h mmH sH h]2mH sH hQhmH sH 2c%9+:+..K1L144::aBbBEE]G^GyGzGKKKKVNXNgd)*$gd)*$a$gd(.Ygd(.YgdMiEgdQgd0[L1|11111e2222222223,3F3Y3333339444455655^666 7>7?7ļ̩̩̞~vnvf^h gmH sH h mH sH himH sH hmH sH hk\7mH sH hMmH sH hB]mH sH h4|mH sH hM^BhumH sH hM^Bh 5mH sH hkcmH sH h)mH sH hmH sH hM^BhmH sH hM^Bh$KmH sH hjumH sH hM^Bh>mH sH hM^Bh8mH sH $?777778(81888999Y9Z9m9r99999999999!:C::::::::::::;;;(;*;/;J;ظظиذ}uhsmH sH h?dmH sH h% mH sH hmH sH hm>5mH sH hM^BhnmH sH hMmH sH hmH sH h1pmH sH h+mH sH h=ZmH sH h.HMmH sH h;%mH sH h gmH sH hvvtmH sH hE#mH sH hB]mH sH .J;L;<<2<x<<<=`=t=====>>#>^>>>>>>>>M?\?x?|???????\@]@^@@@ļıġġyyqyiyh)mH sH hpHmH sH h amH sH h&$mH sH h-mH sH h~mH sH h!mH sH h% mH sH hVmH sH hM^BhnmH sH hJ(mH sH h"mH sH hmH sH hsmH sH hTyhTy6OJQJmH sH hTyhTyH*mH sH hTymH sH (@@@@@@@@@8AkAlAqAAAAAAAA0BBB`BaBbBBBBBBCCmCoCCCDDDDOExE|EEEеz hb 5 h(.Y5 hqG5 h5J5 hpH5 hy5 hDw5 h5h4|mH sH hM^BhnmH sH h0[h5JmH sH hmH sH h5JmH sH hj,mH sH h mH sH h,4mH sH h&$mH sH h)mH sH h HmH sH ,EEJFiF]G^GaGeGyGzGGGGHHHHHHHHHHHHEIFI]IʾpaN$jh ^hj)EHU]mH sH j*bJ hj)CJUVaJ$jhHho,!EHU]mH sH j6 J ho,!CJUVaJjho,!U]mH sH ho,!]mH sH hM^Bho,!]mH sH hI5mH sH hM^BhI5mH sH h][5mH sH hc'5mH sH hM^Bh~dmH sH hpHmH sH hmH sH h]h5]I^I_I`IgIhIIIIIIIIIIIII KKKVKWKKKKKKKKKKKƷƕymmmmmmaUhM^Bho,!5mH sH ho,!ho,!]mH sH hM^Bho,!]mH sH h5I]mH sH $j( hHho,!EHU]mH sH jJ ho,!CJUVaJ$jQ hHho,!EHU]mH sH jJ ho,!CJUVaJho,!]mH sH jho,!U]mH sH $jhHhj)EHU]mH sH jbJ hj)CJUVaJ KKKKKKKKKALZLdLiLLLLL-MFMVMXMgMiMMMMMMMM÷~u~i~`uiTi~Khh]mH sH hM^Bhz]mH sH hkc]mH sH hM^Bh3w]mH sH h6]mH sH hM^Bha]mH sH hM^Bhc]mH sH h;]mH sH hM^Bh<1]mH sH hM^Bh5}]mH sH hM^Bh(5mH sH h5hI5mH sH h5h55mH sH h5hY85mH sH h5h;5mH sH h5h][5mH sH M N`NbN~NNNNNN*O,O.OHOROTOVOXOtOvOzOOOOOOOP"P$P&P(PQQ4Q6Q8Q:QQQRRRRRнвۤММнДۉ{н{{hM^Bh?b6H*mH sH hg#6H*mH sH hg#mH sH hOmH sH hM^Bh$6H*mH sH hM^Bh$mH sH hM^Bh9qmH sH hhmH sH hM^Bh?bmH sH hM^Bh?b6mH sH hM^Bh?b]mH sH hM^Bh@]mH sH ,RRRRRRS S0S2S4S6S>T@TDT,U.U0U8U:UU$V&V^V`VfVVVWW WƸѓxxpepp툓WhM^BhoD6H*mH sH hM^Bh]FmH sH hamH sH h)*mH sH hoomH sH hM^BhoDmH sH hM^BhoD6mH sH hM^Bh)*5]mH sH hoD5]mH sH hM^BhoD5]mH sH hC5]mH sH hM^BhoD]mH sH h5mH sH h/bmH sH hOmH sH hM^Bh$mH sH XNRS2S4S\\\\^^^^ldndff.f/f$p7$8$H$^p`a$gdll $7$8$H$a$gdYr 7$8$H$gd7 $7$8$H$a$gdoD$7$8$H$`a$gdoD$a$gdoD $7$8$H$a$gd?b 7$8$H$gd)* WWWWW W&WX X8XPXnXzXXXXXXYYFYpYzYYYY0Z2Z4ZRZTZB[[[[[\\\\\ܾܶܠܘܾܶ܍܅vjhM^BhoDUmH sH hVmH sH hM^BhFmH sH hmH sH h)*mH sH hM^BhoD6H*mH sH hamH sH hOmH sH he`mH sH hM^BhoD6H*mH sH hM^BhoDmH sH hM^BhoD6mH sH hM^BhoDH*mH sH (\\\\\\\\\\\\]]^^^^6^B^d^f^x^^^^^^Ƚ|nf^нOj*(J hK6CJUVaJh lmH sH hamH sH hM^BhoD6H*mH sH hM^BhoD6H*mH sH hM^BhoD6mH sH hqo3mH sH hqo3hA6mH sH hAmH sH hM^Bh9^mH sH hM^BhoDmH sH hVmH sH jhM^BhoDUmH sH !jhVh}EHUmH sH j"J h}CJUVaJ^^^^^^^^^^^^^^^^^^^^^^:_>_`_b_d_j_~_Ĺxl^SKSKSh$zmH sH hM^BhcdmH sH hK6hYr6