|Home | About | Journals | Submit | Contact Us | Français|
Despite major scientific, medical and technological advances over the last few decades, a cure for cancer remains elusive. The disease initiation is complex, and including initiation and avascular growth, onset of hypoxia and acidosis due to accumulation of cells beyond normal physiological conditions, inducement of angiogenesis from the surrounding vasculature, tumour vascularization and further growth, and invasion of surrounding tissue and metastasis. Although the focus historically has been to study these events through experimental and clinical observations, mathematical modelling and simulation that enable analysis at multiple time and spatial scales have also complemented these efforts. Here, we provide an overview of this multiscale modelling focusing on the growth phase of tumours and bypassing the initial stage of tumourigenesis. While we briefly review discrete modelling, our focus is on the continuum approach. We limit the scope further by considering models of tumour progression that do not distinguish tumour cells by their age. We also do not consider immune system interactions nor do we describe models of therapy. We do discuss hybrid-modelling frameworks, where the tumour tissue is modelled using both discrete (cell-scale) and continuum (tumour-scale) elements, thus connecting the micrometre to the centimetre tumour scale. We review recent examples that incorporate experimental data into model parameters. We show that recent mathematical modelling predicts that transport limitations of cell nutrients, oxygen and growth factors may result in cell death that leads to morphological instability, providing a mechanism for invasion via tumour fingering and fragmentation. These conditions induce selection pressure for cell survivability, and may lead to additional genetic mutations. Mathematical modelling further shows that parameters that control the tumour mass shape also control its ability to invade. Thus, tumour morphology may serve as a predictor of invasiveness and treatment prognosis.
In a healthy body, cells control their proliferation and programmed cell death (apoptosis) in the various tissues so as to optimize body repair and healing. In cancer, this carefully regulated mechanism breaks down; cells proliferate or refrain from dying, may change the microenvironment to favour their survival and may migrate and metastasize in regions far from the primary tumour. This process eventually may kill the host body due to physical obstruction or organ malfunction. Normal cell behaviour is carefully orchestrated through expression of genes within cells and regulatory networks; in cancer, genes that promote proliferation (oncogenes) as well as apoptosis (tumour suppressor genes) may malfunction, and regulatory signals may be ignored. Within an abnormal cell population, additional mutations and epigenetic changes may further lead to different subgroups of cells (‘clones’) that differ in their characteristics. As cells accumulate to form microscopic nodules without access to the vascular network, they receive nutrients and growth factors via diffusion through the neighbouring host (healthy) tissue. Consequently, these nodules typically remain small and grow at most to a few millimetres in diameter.
The accumulation of tumour cells may cause acute and chronic lack of oxygen (leading to hypoxia) and nutrient (e.g., glucose, leading to hypoglycaemia) as well as accretion of metabolites (e.g., lactic acid, leading to acidosis) [226, 257, 265, 532, 550]. As tumour cells accumulate, the insufficiency of the existing vasculature to deliver oxygen and nutrients to all the cells present may induce neovascularization. Cells under stressful conditions will release pro-angiogenic growth factors to drive angiogenesis—the process by which existing blood vessels are stimulated to grow from the main circulatory system to feed tissue with blood, similar to what normally occurs during wound healing. This provides the tumour with a direct supply of nutrients and growth-promoting factors. Once a tumour is vascularized, it can grow larger and even shed cells into the vessels, leading to satellite tumours in distant parts of the body (metastases). Metastasis is the predominant cause of mortality due to cancer. By the time a tumour reaches a clinically detectable size, it is usually in the vascular growth phase. Thus, the transition to metastasis and malignancy typically starts with angiogenesis.
Hypoxia, hypoglycaemia and acidosis are exacerbated by the tumour-induced microvasculature, which, unlike the normal wound healing vasculature tends to be highly disorganized and poorly functioning [258, 291], resulting in considerable heterogeneity in oxygen and nutrient delivery and metabolite removal [289, 290]. These conditions correlate with poor clinical outcome and increased risk of cancer spread through the body [75, 76, 263, 264, 504], may select for apoptosis-resistant tumour cells , induce further blood vessel formation [185, 257, 481], and increase invasiveness [99, 151, 226, 282, 423, 454, 455, 567, 568].
Cell–cell adhesion and communication enable collective, contractile motion by large, multicellular aggregates that move as a functional unit. Collective cell migration dominates in tumours of highly differentiated tissues (e.g. lobular breast cancer and epithelial prostate cancer [56, 72, 83, 401], where invasion by individual cancer cells is rarely detected [198, 425]. It is generally observed  that pre-metastatic phenotypic transitions (e.g. epithelial–mesenchymal transitions, which turn an attached epithelial cell into a motile mesenchymal cell) follow a collective-migration stage and are regulated by the environment (e.g. local hypoxia) [483, 529]. Collective migration has numerous advantages, including high autocrine pro-migratory and protease concentrations and protection of inner cells from immunological defenses. In addition, highly migratory cells at the periphery can promote invasion by less motile but more death-resistant clones. This organized, structured motion of a whole cluster of cells that emerges from the forces and mechanisms that regulate individual cell motion and synchronization suggests that there is more to a tumour than the sum of its parts. In pre-metastatic, locally invasive tumours, collective cell migration in protruding sheets, strands and detached clusters (e.g., pre-metastatic, locally-invasive epithelial tumours) relies mostly on cell–cell and cell–matrix adhesion, cell–cell communications, and the controlled release of matrix-degrading enzymes [196, 198, 425, 461]. The latter degrade proteins such as those that make up the basement membrane separating epithelial cells from stromal tissues, facilitating tumour cell invasion.
The complex interplay of molecular- and tissue-scale dynamics in a tumour can have profound, unforeseen effects on invasion and outcome of therapy. Invasion introduces barriers to resection and treatment; even when resection is an option, it may be challenging to define optimal mass removal because of the difficulty in obtaining an accurate estimate of tumour margin and invasion potential. While angiogenesis plays a critical role in promoting tumour growth and invasion, anti-angiogenic therapy does not always increase length of survival [62, 68, 322, 460], in spite of clinically observed tumour regression [281, 453]. Anti-angiogenic therapy may also exacerbate hypoxia  and cause tumour mass fragmentation, cancer cell migration and tissue invasion [57, 323, 327, 401, 411, 459, 474]. However, pharmacological promoters of cell adhesion are being employed in anti-invasive therapy [25, 52, 153, 168, 261, 278, 326, 344, 370, 566] with inconsistent results , perhaps partly due to drug-induced cellular plasticity [198, 556].
Although molecular mechanisms and cell-scale migration dynamics are well described, the variable empirical and qualitative observations of tumour invasion and response to therapy illustrate the critical need for biologically realistic and predictive multiscale mathematical models that integrate tumour proliferation and invasion with microvascular effects and microenvironmental substrate gradients. Such complex systems, dominated by large numbers of processes and highly nonlinear dynamics, are difficult to approach by experimental methods alone and can typically be better understood with appropriate mathematical models and sophisticated computer simulations, in addition and complementary to experimental investigations. The complexity of cancer progression necessitates a 3D multiscale modelling approach to produce predictive tumour simulators that couple processes occurring at various length and time scales to capture the realistic dynamics involved in the biological environment. While there are more than 100 types of cancer, each with many subtypes, it has been hypothesized that nearly all cancers develop a common set of basic characteristics ; (1) self-sufficiency in growth signals, (2) insensitivity to anti-growth signals, (3) evasion of apoptosis, (4) limitless replicative potential, (5) sustained angiogenesis and (6) tissue invasion and metastasis. By focusing on these common elements, mathematical modelling aims to contribute to the prevention, diagnosis and treatment of this complex disease.
The modelling has the potential to provide important insight into the root causes of solid tumour invasion and metastasis and aid in the understanding of experimental and clinical observations, which may sometimes seem contradictory, and in the design of new, targeted experiments. Biophysically justified mathematical models may also be useful in assessing various treatment strategies. The ultimate goal is for modelling and simulation to aid in the development of individualized therapy protocols to minimize patient suffering while maximizing treatment effectiveness. In order to achieve this goal, mathematical and computational models need to quantify the links of 3D tumour tissue architecture with growth, invasion and the underlying cellular/microenvironmental characteristics.
Numerous mathematical models have been proposed to study the various phases of cancer progression (e.g. see the reviews Adam , Chaplain , Bellomo and Preziosi (2000), Moreira and Deutsch , Bellomo et al , Swanson et al , Araujo and McElwain , Mantzaris et al , Friedman , Ribba et al , Quaranta et al , Hatzikirou et al , Nagy , Byrne et al , Fasano et al , Jackson et al , van Leeuwen et al , Roose et al , Graziano and Preziosi , Harpold et al , Friedman et al , Sanga et al , Martins et al , Deisboeck et al , Anderson and Quaranta , Bellomo et al , Cristini et al , Wang and Deisboeck , Preziosi and Tosin , Jackson et al  and Tracqui ). Most models fall into two broad categories, based on how the tumour tissue is represented: discrete cell-based models and continuum models. Although continuum and discrete approaches have each provided important insight into cancer-related processes occurring at particular length and time scales, the complexity of cancer and the interactions among the cell- and tissue-level scales ideally call for a multiscale continuum-discrete (hybrid) approach, coupling biological phenomena from the molecular and cellular scales to the tumour scale (e.g. see the recent work by Kim et al , Stolarska et al  and Bearer et al ).
In discrete modelling, individual cells are tracked and updated according to a specific set of biophysical rules. This approach is particularly useful for studying carcinogenesis, natural selection, genetic instability and interactions of individual cells with each other and the microenvironment. Analyses of cell population dynamics have also been employed in order to study biological characteristics applying to all cells in a particular population, such as response to therapy and in studies of immunology. There are two main types of discrete models, lattice-based and lattice-free. The former describes the dynamics of discrete tumour cells as automata on a grid whose states are governed by a set of deterministic or probabilistic rules. The latter type describes the actions of discrete cells in arbitrary locations and their interactions. Because these methods are based on a series of rules for each cell, it is possible to translate detailed biological processes (e.g. cell-cycle events, mutation pathways) into rules for the model. On the other hand, the computational cost increases rapidly with the number of cells modelled, limiting these methods in the spatial and temporal scales they can achieve. As a result, a full simulation of a 1 mm3 tumour spheroid, which may contain several hundred thousand cells, is not currently feasible using a solely discrete approach. Further, while these models are capable of describing biophysical processes in significant detail, it may be nontrivial to obtain reliable measurements of model parameters through experiments that can measure the necessary detail at the cell scale.
In larger scale systems, continuum methods provide a good modelling alternative. Continuum models treat tumours as a collection of tissue, where densities or volume fractions of cells are described. This approach draws upon principles from continuum mechanics to describe cancer-related variables as continuous fields by means of partial differential and integro-differential equations. Individual cells and other elements are not tracked; model variables may include cell volume fractions and cell substrate concentrations, e.g. nutrient, oxygen and growth factors. In addition, fast numerical solvers can be developed. Continuum model parameters may be somewhat easier to obtain, analyse and control compared with the discrete case. They may also be more accessible through laboratory experimentation. Although these models are appropriate at the tissue scale where gross tumour behaviour can be quantified, the limitations in scale prevent them from modelling individual cells and discrete events (e.g. epithelial to mesenchymal phenotypic transition that leads to individual cell migration). This may be important when studying the effect of genetic, cellular and microenvironment characteristics on overall tumour behaviour.
Continuum–discrete models attempt to combine both continuum and discrete descriptions of cancer biology in order to bridge the subcellular- and cellular-scales to the tumour scale. In one approach, which we term a composite approach, substances such as oxygen, nutrient, drug, growth factors and certain tissue features (e.g. extracellular matrix) may be described as continuum fields in the tumour microenvironment, while individual discrete elements (e.g. cells or parts of cells) evolve dynamically in response to local conditions such as substance concentration. This modelling strategy suffers from the same limitations as mentioned above for discrete modelling. In a more promising approach, the tumour tissue itself is modelled containing both discrete (cell-scale) and continuum (tumour-scale) elements. This method, which we term a hybrid approach, has the potential to combine the best features of both discrete and continuum models, although more work is necessary to make it competitive with the continuum results obtained at large scales.
In this paper, we provide a limited overview of the theoretical modelling of cancer. While we briefly review discrete modelling, our focus is on the continuum approach, and we discuss hybrid-modelling frameworks where the tumour tissue is modelled containing both discrete (cell-scale) and continuum (tumour-scale) elements, thus connecting the micrometre to the centimetre tumour scale. We also review recent examples that incorporate experimental data into model parameters. We limit the scope further by considering models of tumour progression that do not distinguish tumour cells by their age. We do not consider tumour immune system interactions nor do we describe models of therapy. For age-structured models, see, for example, [2, 46, 165, 166, 537, 538], and for tumour immune system modelling see, for example, [6, 7, 58, 61, 147, 184, 354, 516]. Further details and references may be found in the review papers listed above.
As we show, mathematical modelling predicts that transport limitations of cell nutrients, oxygen and growth factors may result in cell death that leads to morphological instability, providing a mechanism for invasion via tumour fingering and fragmentation. These conditions induce selection pressure for cell survivability, and may lead to additional genetic and phenotypic changes. Mathematical modelling shows that parameters that control the tumour mass shape also control its ability to invade. Thus, tumour morphology may serve as a predictor of invasiveness and treatment prognosis.
The outline of this paper is as follows. In section 2 we review continuum modelling and the incorporation of biologically relevant parameter values into multiscale models of tumour growth and invasion. A basic model founded on classical work is reviewed first, and then expanded to tumour growth in heterogeneous tissue and the host vascular response (angiogenesis). We then consider multiphase modelling to simulate multiple cell species, and include effects from cell adhesion and chemotaxis. In particular, we evaluate the relevance of theoretical cancer modelling to patients suffering from glioma (brain) tumours and ductal carcinoma in situ (DCIS) (breast) tumours. We briefly review discrete modelling in section 3. In section 4 we discuss hybrid modelling, where the tumour is described using both continuum and discrete elements, and which is capable of connecting cell and tissue scales to provide practical as well as theoretical insight into cancer growth. Conclusions and future directions are described in section 5.
Continuum tumour models are based on reaction–diffusion equations describing the tumour cell density (e.g. [23, 544, 414])), extracellular matrix (ECM), matrix-degrading enzymes (MDEs) (e.g. [75, 76, 111, 264, 287]) and concentrations of cell substrates such as glucose, oxygen, and growth factors and inhibitors (e.g. [69, 104, 109, 225, 257, 263, 353, 391]). Classical work [250, 251] used ordinary differential equations to model tumours as a homogeneous population, as well as partial differential equation models confined to a spherical geometry. In the case of avascular tumours, growth has been modelled as a function of cell substrate concentration, usually oxygen. More recent work has incorporated cell movement, through diffusion (e.g. [80, 109, 130, 293, 478, 480, 508–511, 557]), convection (e.g. [93, 139, 414, 543, 546]) and chemo/haptotaxis (e.g. [364, 414, 478]). Cell proliferation, death and pressure have also been considered (e.g. [3, 22, 39, 59, 73, 84–87, 92, 93, 95, 96, 120, 129, 146, 186, 189, 220, 222, 272, 283–285, 287, 298, 328, 336, 346, 414, 420, 457, 468, 493, 496, 521, 536, 545]). Linear and weakly nonlinear analyses have been performed to assess the stability of spherical tumours to asymmetric perturbations (e.g. [8, 34, 88, 87, 89, 91, 92, 96, 109, 129, 235, 339]) in order to characterize the degree of aggression. Various interactions with the microenvironment, such as nutrient-, inhibitory factor- or stress-induced limitations to growth, have also been studied (e.g. [4, 5, 19, 21, 23, 35, 37, 38, 107, 108, 138, 298, 457]). The models may account for observations of stronger cell–cell interactions (cell–cell adhesion and communications), high polarity and strong pulling forces exchanged by cells and ECM ([117–119, 198]). ECM reorganization by tumour cells  has been incorporated, and various degrees of dependence of the cells on signals from the matrix have been modelled. The models are typically single species (e.g. single-phase tumours), treating the tumour (or more generally biological tissues) as fluid (e.g. [73, 74, 89, 91, 92, 95, 105, 206, 251]), elastic/hyperelastic (e.g. ( [18, 21, 35, 210, 211, 218, 279, 298, 372, 477, 536]), poroelastic (e.g. ), viscoelastic (e.g. [346, 317], and elasto-viscoplastic (e.g. . Theoretical nonlinear analyses of various tumour models mentioned above have been performed (e.g. [53, 78, 122, 131–134, 136–142, 155, 182, 199, 202–204, 207, 208, 389, 404, 480, 515, 517, 539, 558, 562, 564, 576, 577]). Recently, mixture models have been developed that are capable of describing the detailed interactions among multiple solid (cell) species and extra-/intra- cellular liquids (see section 2.5).
In the next subsections, we first describe tumour growth in homogeneous tissues by posing a basic continuum tumour model, focusing on the formulation of [129, 339] based on the classical work mentioned above. Model parameters are calibrated from experimental data. Later, in section 2.3, we discuss tumour growth in complex, heterogeneous tissues. Tumour mechanical responses are also reviewed. The basic model is then extended in section 2.4 to include angiogenesis and vascular growth.
A basic tumour model represents tumour cells growing as a sphere-like structure without direct access to the vasculature. During this avascular growth, tumour cells receive oxygen, nutrients and growth factors via diffusion through the host tissue. This phase can be investigated by in vitro experiments where cancer cells are cultured in a three-dimensional geometry [315, 324, 383, 384, 505, 506, 540]. Due to cell–cell adhesion, certain cancer cell lines will self-organize into multicellular, roughly spheroidal colonies. The outer cells tend to proliferate while the cells in the interior necrose (die) due to lack of nutrients. For example, the typical distance an oxygen molecule will diffuse before being uptaken is approximately 100 μm. This limits the size to which a tumour spheroid can grow (1–2 mm in diameter). A layer of quiescent (hypoxic) cells separates the necrotic core from the proliferating rim. Because of the three-dimensionality, the growth of multicellular spheroids may be similar to that of in vivo avascular tumours. There is a significant amount of experimental data in the literature on the internal structure of multicell spheroids and the spatio-temporal distribution of cell substrates (see references above). Thus, this is a good model system to test mathematical predictions.
Greenspan [250, 251] developed one of the earliest continuum models of tumour growth as a function of diffusion of cell substrates, as observed in previous studies (e.g. [81, 518]). Shymko and Glass  accounted for a mitotic inhibitor and analysed the stability of growth, McElwain and Morris  accounted for apoptosis, and Adam  discussed the immune response. Byrne and Chaplain [89, 91, 93] studied the growth and stability of radially symmetric tumours without and with necrosis, as well as the effects of cell substrates and inhibitors. Chaplain  presented mathematical models of spherical tumour growth through the stages of avascular growth, angiogenesis and vascularization, as well as pattern formation in cancer . Friedman and Reitich [206, 207] and Cui and Friedman  studied nonnecrotic vascularized radially symmetric and spatially patterned tumours modelled through a free-boundary problem, where tumour growth was dependent on the level of diffusing cell substrates. See also [131, 134, 136, 202, 253, 515, 517, 558] for later extensions to a variety of different tumour models.
Based on these and other classical continuum tumour models, Cristini et al  performed computer simulations of tumour growth beyond the limited capabilities of mathematical linear analyses and spherical geometries, thus enabling the nonlinear modelling of complex tumour morphologies. Using a new formulation of these classical models, they showed that tumour evolution could be described by a reduced set of two dimensionless parameters (related to mitosis rate, apoptosis rate, cell mobility and cell adhesion), independent of the number of spatial dimensions. These parameters regulate the morphology and growth (invasiveness) of avascular and vascularized tumours. Critical conditions were predicted that separate compact, noninvasive mass growth from unstable, fingering, infiltrative progression  thus suggesting that the mechanisms that control tumour morphology also control its ability to invade. Recently, Li et al  have extended this work to arbitrary geometries in 3D. This morphological instability provides a mechanism for invasion without angiogenesis and may allow the tumour to overcome diffusional limitations to growth by creating excess surface area that exposes more interior cells to oxygen and nutrient. Indeed, numerical simulations show that the tumour grows unbounded by repeated sub-spheroid growth (budding), fingering and folding to create a complex shape. That is, morphological instability driven by microenvironmental gradients of nutrient and oxygen selects for locally higher cell proliferation. Tumours may thus escape diffusion-limited constraints without recourse to angiogenesis, as has been observed experimentally (e.g. [148, 194]).
Define Ω(t) to be the tumour domain, Σ(t) to be the boundary between the tumour tissue and the host tissue, n to be the unit outward normal vector to Σ and x to be the position in space. See figure 1. Following Greenspan , Byrne and Chaplain , Friedman and Reitich , Cristini et al , and others, it may be assumed that in the proliferating tumour domain, the cell density is constant. Therefore, mass changes correspond to volume changes. Defining u to be the cell velocity, the local rate of volume change · u is given by
where λp is the cell-proliferation rate and λp is given by
where n denotes the concentration of a cell substrate (e.g. oxygen or glucose). The first term in equation (2) corresponds to the rate of volume growth due to mitosis while the second is the rate of volume loss due to apoptosis (programmed cell death). Here, λA is the rate of apoptosis (which may actually depend on n) and b is a measure of mitosis.
Cell substrates diffuse through the extracellular matrix (ECM), as well as intracellularly, and are uptaken by tumour cells. Since the rate of diffusion of oxygen (or glucose) is much faster (e.g. ~1 min−1) than the rate of cell proliferation (e.g. ~1 day−1), the substrate may be regarded to be in a steady state for a given tumour morphology (e.g. [89, 129, 206, 251]). This gives:
where Γ is the rate at which nutrient is added to Ω and is given by
where the first term describes the source of nutrient from the vasculature while the second describes nutrient uptake by cells. Here, λB is the blood-tissue transfer rate of nutrient, nB is the concentration of nutrient in the blood and λ is the rate of consumption of nutrient by the tumour cells. In this simplified model, the vasculature is assumed to be uniform, and thus vascular growth is associated with a bulk source of oxygen, nutrients and growth factors. Note this implies that growth is limited by the diffusion of the cell substrates.
where P is the oncotic (solid) pressure and μ is a mobility that reflects the combined effects of cell–cell and cell–matrix adhesion. Alternatively, the velocity may be determined using the Stokes equations (e.g. [186, 189, 204, 558]) or the Darcy-Stokes (Brinkman) equations ((e.g. ). Models of viscoelasticity (e.g. ), elasto-viscoplasticity  and soft tissue elasticity may also be employed, see Internal stress in section 2.3.
The boundary conditions on the tumour interface Σ may be taken to be
where the pressure boundary condition (7) reflects the influence of cell–cell adhesion through the parameter γ and κ is the local total curvature. For simplicity, here assume n∞ is constant so that outside the tumour, the nutrient is uniform. Nutrient inhomogeneity in the tumour microenvironment in 2D has been considered in [126, 194, 350, 351, 575] and more recently in 3D in [193, 464, 555].
The normal velocity V = n · (u)Σ of the tumour boundary is
Following [89, 129, 206, 250] and others, assume that λ, λA, λB, nB, b are uniform. Following , denote λM = bn∞ to be the characteristic mitosis rate, to be the intrinsic relaxation time scale, and B = nBλB/n∞(λB + λ) to be a measure of the extent of vascularization. Introducing the non-dimensional length scale , and time scale , a modified concentration and pressure can be defined :
where G and A measure the relative strength of cell–cell and cell–matrix adhesion and apoptosis, respectively:
Dropping the bar notation, the non-dimensional equations for Γ and p can be obtained:
with boundary conditions:
in a d-dimensional tumour (d = 2, 3). The non-dimensional normal velocity of the tumour–host interface is
A study of spherically symmetric tumour growth provides insight into the regimes of growth described by the model (e.g. [89, 129, 339]). In this case, the PDEs reduce to ODEs in the polar coordinate r. Accordingly, from equation (15) the evolution equation for the tumour radius R is
For a radially symmetric tumour, |G rescales time. In all dimensions, unbounded growth (R → ∞) occurs if and only if AG ≤ 0. Three regimes of growth are identified, and the behaviour is qualitatively unaffected by the number of spatial dimensions d.
G ≥ 0 and A > 0 (i.e. B < λA/λM). Note that the special case of avascular growth (B = 0) belongs to this regime. The evolution is monotonic and always leads to a stationary state R∞ (if A > 1, then R∞ = 0). This behaviour is in agreement with experimental observations of in vitro diffusional growth  of avascular spheroids to a dormant steady state [383, 505]. In the experiments, however, tumours always develop a necrotic core that further stabilizes their growth .
G ≥ 0 and A ≤ 0 (i.e. 1 > B ≥ λA/λM). Unbounded growth occurs from any initial radius R0 > 0. The growth tends to exponential for A < 0 with velocity V → −AGR/d as R → ∞, and to linear for A = 0 with velocity V → G as R → ∞.
G < 0 (i.e. B > 1). For A > 0, growth (V > 0) may occur, depending on the initial radius, and is always unbounded; for A < 0 (for which cell apoptosis is dominant: λA/λM > B), the evolution is always to the only stationary solution R∞ = 0. This stationary solution may also be achieved for A > 0. The stationary radius R∞ is independent of G, and is a solution of V = 0 with V from (16).
Consider a perturbation of the radially symmetric tumour interface Σ:
where rΣ is the radius of the perturbed tumour–host interface, δ is dimensionless perturbation size and Yl,m is a spherical harmonic, where l and θ are polar wavenumber and angle and m and ϕ are azimuthal wavenumber and angle, respectively. By solving the system of (11)–(15) in the presence of a perturbed interface and matching powers of the perturbation δ, we obtain the evolution (16) for the unperturbed radius R and the following equation for the shape factor δ/R [129, 339]:
Note that δ/R is the appropriate way to measure the perturbation since the underlying radius of the symmetric tumour is time-dependent. Also observe that the linear evolution of the perturbation is independent of the azimuthal wavenumber m and that there is a critical mode lc such that perturbations grow for l < lc and decay for l > lc. The critical mode depends on the parameters A, G and the evolving radius R. Prior to the work of Cristini et al and Li et al linear analyses [89, 93, 251] considered only the special case where the unperturbed configuration is stationary (i.e. R constant in time).
In the low-vascularization regime, the existence of non-symmetric, steady-state tumour shapes may be predicted by linear theory. This is seen as follows. The stationary radius R∞ is solution of (16) with V = 0. An analysis  of equation (18) reveals that there exists a non-negative critical value G = Gl(R∞, A) such that perturbation also remains stationary. The perturbation δ/R∞ grows unbounded for G > Gl and decays to zero for G < Gl.
The analysis described above for stationary states may be extended [129, 339] to the case in which the underlying symmetric tumour is time-dependent. In particular, one may obtain self-similar growth for a single mode l. This can be done as follows. Take G to be constant and take A = A(l,G,R) such that (d/dt)(δ/R) = 0. This gives,
In figure 2 (top) the apoptosis parameter A(l,G,R) is shown for d = 2 (dashed) and d = 3 (solid) from [129, 339]. The growth velocity V corresponding to=self-similar evolution, obtained from (16) with A = A(l,G,R), is plotted in figure 2 (bottom). The A-curves separate parameter space into regions of stable and unstable growth of a given mode l. This figure also indicates that in the low-vascularization (diffusion-dominated) regime where G > 0 and A > 0), self-similar evolution towards a stationary state is not possible for G constant. Instead, one may obtain self-similarly growing and shrinking tumours. In the moderate and high-vascularization regimes, the only self-similarly evolving tumours shrink to a point.
This work shows that if one is able to control the apoptosis parameter in a subtle way via equation (19) by applying appropriate therapeutic treatments, for example, one has the possibility to prevent the tumour from becoming unstable and invasive. This has important implications for the angiogenic response of the host—a smaller surface area to volume ratio means less flux of angiogenic factors—as well as for the resectability of the tumour—a compact shape is easier to remove surgically.
The linear stability analysis [129, 339] shows that during growth, perturbations can increase only in the low-vascularization regime. In the moderate and high-vascularization regimes, perturbations always decay during growth. For example, taking A to be a non-negative constant, observe from figure 2(top) that the evolution will cross the A(l,G,R) curves and hence becomes unstable only for G > 0. Instability arises because growth is limited by diffusion of nutrient (e.g. diffusional instability). This is analogous to the Mullins–Sekerka diffusional instability that occurs in crystal growth . In the high-vascularization regime, shrinkage may be unstable.
The parameters A and G in the basic tumour model can be estimated by comparison with experimental tumour spheroid results. An estimation was performed by Frieboes et al , using ACBT (grade IV human glioma multiforme) tumour spheroids, where experiments were performed and the results were analysed using the basic mathematical model. In their in vitro experiments, the glucose and fetal bovine serum (FBS) concentrations were varied to modify cell adhesion and rates of cell proliferation. Rates of cell proliferation increased with serum concentration as expected. Cases with 1% FBS had slowest tumour mass growth, whereas 10% FBS cases had the fastest. Note that serum may also alter adhesion while increasing proliferation. Furthermore, cell mobility (adhesion) was found to increase (decrease) with glucose concentration, in agreement with previous observations that higher levels of glucose may reduce oxygenation in large spheroids . Thus, an important effect of higher glucose was to decrease cell adhesion forces and thereby contribute to morphologic instability. Most stable, compact morphologies were observed with low/medium levels of both glucose and FBS, in which tumours shed fewest cells and attained smallest overall sizes. In contrast, the combination of high glucose and any serum concentration exhibited very unstable morphologies because cells were very motile. Similarly, the combination of any glucose and 10% FBS had very unstable morphologies apparently because cells proliferated faster than they had time to connect into a stable structure.
By estimating the size of the viable rim of cells on the tumour periphery, the diffusion length LD ≈ 100 μm is obtained. By fitting equation (16) to the spheroid growth curves at small times (nearly exponential growth), the proliferation rate is estimated to be λM ≈ 1 day−1. The calibrated model is consistent with other measurements. For instance, equation (3) gives an oxygen penetration length scale Loxy = (Doxy/λoxy,uptake)−1/2. By measuring the distance between the necrotic core and the basement membrane, this length can be estimated to be 100 to 140 μm; using previously published values λoxy,uptake = 9.41 × 10−2 s−1  and Doxy = 1.45 × 10−5 cm2 s−1  gives Loxy = 124 μm, in excellent agreement with these results. Similar calculations were consistent for calculating the glucose penetration length and uptake rate [230, 252, 301], confirming that hypoxia is the limiting factor for tumour cell viability.
An analytical relation between A and R can be established for a steady-state solution by taking V = 0 in equation (16). The corresponding steady-state relation A = As(R) can be used to estimate the value of A in the experiment by taking R to be the average experimental tumour radius, non-dimensionalized by LD, of morphologically stable spheroids (low/medium levels of both glucose and FBS). The parameter A was thus estimated as 0.26 ≤ A ≤ 0.38. This is an over-estimate in the proliferating rim, however, because the mathematical model described assumes that cell death occurs uniformly throughout the spheroid. In fact, cell-death should be spatially heterogeneous with the largest values occurring in the interior hypoxic region where cells are starved of nutrients.
Frieboes et al (2006) showed that the parameter G can be estimated by comparing the pressure in the proliferating rim with the pressure at the tumour boundary. In the proliferating rim at steady state, the dimensional pressure is approximately while at the tumour boundary P ≈ 2τ/LDR, where R is the non-dimensional (stable) tumour spheroid radius which follows from equation (7). Equating the two and using the definition of G with B = 0, again with R as the average experimental tumour radius, non-dimensionalized by LD, the estimate 0.6 ≤ G ≤ 0.9 is obtained for stable tumour spheroids. Spheroids with values of G above this range will be morphologically unstable due to weak adhesive forces. Interestingly, this approach provides an indirect method for estimating G without directly measuring cell–cell adhesion .
Given A and G, Frieboes et al (2006) used the linear stability analysis to predict the morphological stability of the tumour spheroids. In figure 3  a phase diagram is presented in which the marginal stability curves A(l,G,R) from equation (19) are plotted with l = 4 and G = 0.9 and G = 0.6. The parameter l = 4 was chosen because unstable morphologies in vitro seemed to exhibit mostly tumour surface perturbations characterized by low wave numbers (e.g. 3 or 4) at the onset of the instability. The stationary relation As(R) is also shown. The horizontal lines indicate the experimentally estimated values A = 0.26 and A = 0.38 obtained as described above. The experimental range of the parameters A and G is indicated by the shaded region. In the presence of cell substrate gradients, morphology can be ‘unstable’ when cell adhesion is weak (large G), whereas for small G, tumour morphology is ‘stabilized’ by cell adhesion . The larger a tumour grows, the weaker the stabilizing effect of cell adhesion. Each A(l,G,R) curve describes a tumour with specific cell characteristics and divides the parameter space into stable (on the left) and unstable regions (on the right). The lower the cell adhesion, the more shifted to the left the G curve is, thus reducing the range of sizes of tumours that will be morphologically stable. As a tumour grows, this corresponds to moving from left to right and thus may lead to eventually crossing the G-curve corresponding to that tumour's degree of morphological stability.
In figure 3, the filled symbols denote experimental spheroids. The spheroid denoted by the filled square is very stable, being on the left of all curves that are compatible with stable morphologies; the spheroid denoted by the filled diamond is very unstable and the spheroid denoted by the filled circle is marginally stable (i.e. it may develop morphologic instability depending on the value of G). In their experiments, Frieboes et al (2006) were able to observe both stable and unstable spheroids by varying the parameter G. This was achieved by altering the glucose and FBS concentrations consistent with the theory since FBS and glucose affect the parameter G by affecting cell proliferation and adhesion as described above. Thinking of this another way, this work shows that by training the model on stable spheroid data to estimate G, the model is able to predict both stable and unstable spheroid morphologies.
To investigate the effect of nonlinearity, nonlinear analyses have been performed. A sample of results is as follows (see also references above). The well posedness of solutions has been proven by Friedman and Reitich  and Cui and Friedman  for radially symmetric solutions using the Darcy law model; Cui  for a model with inhibitors; Friedman and Reitich , Bazaliy and Friedman [53, 54], Cui , Cui and Wei  and Cui  for non-symmetric solutions; Friedman  and Wu and Cui  for a Stokes model; Zhou, Escher and Cui  for a model of multilayer tumours and Cui  for a hyperbolic equation model of tumour growth. Bifurcations from spherically symmetric solutions have been studied by Friedman and Reitich , Friedman and Hu [203–205] and Zhou and Cui . In recent work, Xu and Gilbert  formulated and analysed inverse problems that connect a tumour model for ductal carcinoma in situ with clinical data; the authors also performed two-dimensional numerical simulations.
Because of the difficulty in rigorously analysing the models, and the limitations of theoretical analyses (e.g. proofs are often not constructive), numerical simulation is critical to obtaining insight into the behaviour of solutions in the nonlinear regime. For the Darcy law model presented earlier, efficient numerical algorithms have been developed (e.g. Cristini et al  and Li et al ) to solve equations (11)–(15) in two- and three-dimensions. In the numerical approach, the partial differential equations in the whole domain were reformulated into boundary integral equations that hold only on the tumour/host interface, using potential theory.
Friedman and Reitich [207, 208] proved that there exist non-symmetric steady tumour shapes that are solutions of the fully nonlinear equations. Their proof was not constructive, however. In , the numerical scheme described above was used to obtain approximations of these solutions in 2D. In the nonlinear case, non-symmetric steady tumour shapes may be found by taking , where Gl is the result from linear theory (section 2.2.4). It is found that , where δ is a measure of the perturbation size. Thus, nonlinearity is destabilizing for the stationary shapes.
The effect of nonlinearity on the self-similar evolution for d = 2 predicted from the linear analysis by Cristini et al (2003) can be investigated. As discussed in section 2.2.4, linear self-similar evolution requires the time-dependent apoptosis parameter A = A(l,G,R). In the nonlinear regime, Cristini et al (2003) used this linear relation with the radius R, determined from the area of an equivalent circle: . In figure 4, linear (dashed) and nonlinear (solid) solutions are compared in the low-vascularization regime for l = 5, G = 1, A = A(l,G,R) and R0 = 4. Since V < 0, the tumours shrink and A increases. In the left frame, the initial perturbation is δ0 = 0.2 and in the right the perturbation initially is δ0 = 0.4. Results reveal that large perturbations are nonlinearly unstable and grow, leading to tumour fragmentation. This can have significant implications for therapy. For example, one can imagine an experiment in which a tumour is made to shrink by therapeutic intervention such that A is increased by increasing the apoptosis rate λA. This example shows that a rapid decrease in size can result in shape instability leading to tumour break-up and the formation of microscopic tumour fragments that can enter the blood stream through leaky blood vessels, thus leading to metastases.
All the nonlinear simulations of growth in the high-vascularization regime lead to stable evolution, in agreement with the linear analysis, i.e. well-vascularized tumours tend to grow in compact, nearly spherical shapes showing little or no signs of invasiveness. This prediction suggests that tumours could maintain stable morphology under more normal microenvironmental conditions, as has been observed in experiments [57, 323, 327, 390, 459, 490].
In some cases, it has been experimentally observed that highly vascularized cancers evolve invasively by extending branches into regions of the external tissue where mechanical resistance is lowest (e.g. ). These results suggest that formation of invasive tumours should be due to anisotropies rather than to vascularization alone. Anisotropies (e.g. in the distribution of the resistance of the external tissue to tumour growth, or in the distribution of blood vessels) are neglected in the model studied here but will be included in the next section. This consequence, which had not been recognized before, is supported by recent experiments  of in vivo angiogenesis and tumour growth.
In figure 5, the evolution of the 2D tumour surface from a nonlinear boundary integral simulation (solid) is compared with the result of the linear analysis (dotted), using A = 0.5, G = 20 . According to linear theory (equation (16)), the tumour grows. The radially symmetric equilibrium radius R∞ ≈ 3.32. Mode l = 2 is linearly stable initially, and becomes unstable at R ≈ 2.29. The linear and nonlinear results in figure 5 are indistinguishable up to t = 1, and gradually deviate thereafter. A shape instability develops and forms a neck. At t ≈ 1.9 the linear solution collapses suggesting pinch off. However, the nonlinear solution is stabilized by the cell-to-cell adhesive forces (surface tension) that resist development of high negative curvatures in the neck. This is not captured by the linear analysis. Instead of pinching off, as is predicted by linear evolution, the nonlinear tumour continues to grow and develops large bulbs that eventually reconnect thus trapping healthy tissue (shaded regions in the last frame in figure 5) within the tumour. The frame at t = 2.531 describes the onset of reconnection of the bulbs. It is expected that reconnection would be affected by diffusion of nutrient outside the tumour, which is not included in the model used here (but will be considered in the next section).
An analogous evolution is observed in 3D . See figure 6, where two 3D views of the morphology are shown. The tumour does not change volume in the simulation because of the spatial rescaling. At early times, the perturbation decreases and the tumour becomes sphere-like. As the tumour continues to grow, the perturbation starts to increase around time t ≈ 0.4. The tumour begins to take on a flattened ellipse-like shape. Around time t ≈ 2.2 the perturbation growth accelerates dramatically, and dimples form around time t ≈ 2.42. The dimples deepen, and the tumour surface buckles inwards.
The biological significance of the diffusional instability  is that the instability allows the tumour to increase its surface area relative to its volume, thereby allowing the cells in the tumour bulk greater access to nutrient. This in turn allows the tumour to overcome the diffusional limitations on growth and to increase to larger sizes than would be possible if the tumour were spherical. Thus, diffusional instability provides an additional pathway for tumour invasion that does not require an additional nutrient source such as would be provided from a newly developing vasculature through angiogenesis. These predictions of development of shape instabilities are in agreement with experimental observations (e.g. [196, 198]).
Microenvironmental inhomogeneities play a significant role in the growth of a tumour [171, 263, 411, 466]. For example, hypoxic microenvironments induce tumour and endothelial cells to upregulate HIF-1 target genes leading to the secretion of proangiogenic factors, matrix degrading enzymes and decreased cell–cell and cell–matrix adhesion [178, 306, 424]. In addition, hypoxic microenvironments affect the metabolism of tumour cells leading to activation of the glycolytic pathway and acidosis in the microenvironment [226, 228, 231]. The presence of different tumour and host tissues (e.g. grey/white matter in the normal brain) and variations in the extracellular matrix also influence tumour progression through chemical and physical interactions . In addition, the microenvironment may induce clonal diversity .
A key feature of the hypoxic tumour microenvironment is the existence of oxygen gradients. These gradients may arise from inadequate vascularization, exacerbated by disordered tumour-induced angiogenesis (see section 2.4). This may induce necrosis. The basic model in section 2.2.2 was extended by Byrne and Chaplain  to include necrosis. There have been many studies of the effects of necrosis on tumour growth. See the review papers listed earlier. For example, Garner et al  incorporated necrosis in a model of spherical tumour growth and used the conservation of energy to obtain scaling laws for the growth of the tumour and necrotic core. In [348, 350, 351] a further extension was introduced to study the effects of oxygen variation in the tumour microenvironment. In this model, an avascular tumour is modelled to occupy volume ΩT(t) with boundary Ω, denoted by Σ, viable region ΩV and necrotic region ΩN where tumour cells die due to low nutrient levels. The viable region is divided into a proliferating region ΩP where the nutrient levels are high enough to permit cell proliferation and a quiescent/hypoxic region ΩQ where nutrient levels are insufficient to sustain proliferation. The growing tumour also interacts with the surrounding microenvironment in the host tissue; this region is denoted by ΩH. In addition, equation (3) governing the distribution of oxygen (or any other vital cell substrate) is generalized to allow for non-constant diffusivity and variable uptake and natural decay. Using the same non-dimensional scaling as before, this gives
where λn is the uptake/decay rate. It is assumed that oxygen is uptaken at different rates in the different domains of the tumour and the host tissue (ΩH). For example, the uptake may be different in domains where the cells are proliferating ΩP or quiescent ΩQ. In addition, oxygen may be degraded in the necrotic domain ΩN due to the presence of oxidizing agents . These domains may actually be defined in terms of the oxygen concentration as in equation (21). In equation (20), is the nutrient delivery rate from the pre-existing vasculature, B(x) is the pre-existing vascular density and χH is the characteristic function of the host domain ΩH. Normalized by the nutrient uptake in the proliferating tumour region, the uptake and decay function and definitions of the domains are taken to be
where np and nN are the nutrient thresholds for proliferation and necrosis, respectively. This function reflects the fact that nutrient is uptaken much faster in the tumour than in the host tissue (hence the relative uptake in ΩH is zero) and the fact that when cells necrose, they release their cellular contents which are oxygen reactive [212, 318] and thus this effect on the nutrient concentration can be modelled through the decay rate . Note that with this choice of uptake/decay term, the nutrient equation is nonlinear.
Cells and ECM in the host tissue ΩH and viable tumour region ΩV are affected by a variety of forces, each of which contributes to the cellular velocity field u. Proliferating tumour cells in ΩV generate an internal (oncotic) mechanical pressure (hydrostatic stress) that also exerts force on surrounding non-cancerous tissue in ΩH. Tumour, non-cancerous cells and ECM can respond to pressure variations by overcoming cell–cell and cell–ECM adhesion and moving within the ECM scaffolding of matrix proteins that provides structure to both tumour and host tissue.
Using Darcy's Law as a constitutive relation for the cell velocity in the host and tumour domains, and a domain-dependent net cell-proliferation rate λp, the mechanical pressure is given by
where the right-hand side is λp. Note that the cellular mobility μ also measures permeability of tissue to tumour cells. See [23, 95] for further motivation of this approach from the perspective of mixture modelling. In the host domain, the non-dimensional net proliferation rate is assumed negligible since tumour cells proliferate more rapidly or die at a lower rate than host cells. In the proliferating domain, the proliferation rate is assumed to be linear in the oxygen concentration and apoptosis is assumed to result in volume loss as in section 2.2.2. The parameter GN is the non-dimensional rate of volume loss due to necrosis [348, 575]. Following  and others, cell–cell adhesion forces are modelled in the tumour by generalizing the condition (7) and introducing a Laplace-Young jump condition:
where κ is the mean curvature and . It is also supposed that the pressure jump across the necrotic interface is zero, which reflects the low cell–cell adhesion in the perinecrotic region and the increased cellular mobility observed in hypoxic cells [75, 99, 265, 423, 424, 454]. The cell velocity is also assumed to be continuous across the tumour/host interface and the boundary of the necrotic tumour domain. Cellular proliferation and death are in balance in the far-field, i.e.
In Macklin and Lowengrub , an accurate ghost-cell [183, 244]/level-set  method was developed to solve the system described above. In this approach, the tumour–host interface is described as the zero level-set of an auxiliary function (level-set function). The equations are discretized on a Cartesian mesh and the stencil is adapted by introducing ghost cells near the interface to achieve accurate results. This algorithm is capable of describing complex morphologies evolving in heterogeneous domains. Jumps in the normal derivative are discretized without smearing jumps in the tangential derivative. A new adaptive solver for linear and nonlinear quasi-steady reaction–diffusion problems (NAGSI), an adaptive normal vector discretization for interfaces in close contact, and an accurate discrete approximation to the Heaviside function were introduced.
The effects of the tumour microenvironment on the morphology and growth patterns of 2D can be considered by modelling avascular tumours growing into piecewise homogeneous tissues. Parameters in the model can be varied to represent these effects in a simplified fashion. For example, oxygen gradients in the microenvironment can be created by varying the ratio of oxygen diffusivities in the host domain ΩH and the tumour domain ΩT (the diffusivity in each region is assumed to be uniform). The permeability of the host domain to the tumour cells can be varied through the ratio of respective mobilities (also assumed to be uniform in each domain). The actual situation can be more complex. For instance, when tumour cells are hypoxic, cellular pathways that stimulate cell migration may be activated [178, 263, 306, 331, 424]. This may be modelled by increasing the mobility μ as oxygen level decreases [352, 555] or as a chemotactic response to oxygen gradients , see below and section 2.5.4.
Using the diffusion and mobility parameters as described above, the effects of the tumour microenvironment on growth can be characterized through a morphology diagram as shown in figure 7 . Growth is simulated over a wide range of microenvironmental parameters (with D and μ being the ratios of diffusion and mobility parameters in ΩH and ΩT) with G = 20; GN = 1 and nN = 0.35, each with identical initial shape. The apoptosis parameter A = 0 because on the time scale considered the tumours apoptosis was assumed not to occur. Further, the quiescent region here is not considered, i.e. nP = nN.
In figure 7, the shape of each tumour is plotted at time T = 20.0 (scaled with the mitosis time ~1 day so that this corresponds to approximately 20 days). In all figures, the black regions denote ΩN where the tumour is necrotic, the grey regions show the viable tumour region ΩV, and the white regions correspond to ΩH, which consists of the ECM, non-cancerous cells, and any other material outside of the tumour.
On the horizontal axis, the nutrient diffusivity of the surrounding tissue is varied; as D increases from left to right, the simulated microenvironment varies from nutrient poor (left, hypoxic) to nutrient rich (right, normoxic). On the vertical axis, the mobility of the surrounding material is varied; as μ increases from bottom to top, the microenvironment ranges from low mobility (bottom) to high mobility (top). The greater the mobility μ, the greater the ability of the external, non-cancerous tissue to respond to the pressure generated by the growing tumour and the easier it is for tumour cells to invade the host tissue.
Three distinct tumour morphologies were observed through this range of simulated tissue types. In the nutrient-poor regime on the left side of the diagram, the growing tumours break into fragments so as to optimize the number of tumour cells with access to nutrient. The nutrient-rich, low-mobility regime in the bottom right of the morphology diagram is characterized by invasive fingering, where buds develop on the tumour boundary and elongate. This increases the surface area to volume ratio and enables interior tumour cells to have greater access to nutrient (which is widely available in the microenvironment). The nutrient-rich, high-mobility regime in the top right of the diagram exhibits compact/hollow growth, where the tumours tend to grow into spheroids and typically form abscesses filled with non-cancerous tissue, cellular debris and fluid, similar to a necrotic core. These morphologies are similar to those observed experimentally in vitro (e.g. [194, 350]).
Tumour morphologies in figure 7 are qualitatively similar when recomputed with different genetic/phenotypic characteristics (modelled by G, GN and σN), although large changes in the genetic/phenotypic parameter values can shift the morphology from one type to another. Thus, a tumour's morphology seems to critically depend upon the characteristics of the microenvironment.
Tumours growing into nutrient-poor microenvironments demonstrate repeated fragmentation through a wide range of mitosis rates (governed by the parameter G) and necrotic tissue degradation rates (GN). Tumour fragmentation was observed in the simulations in almost all cases, particularly for fast-proliferating, aggressive tumours with higher values of G. Similarly, increasing the rate of necrotic tissue degradation (GN) tends to destabilize the tumour, also leading to an increased rate of fragmentation. However, this effect is highly nonlinear: if GN is large relative to G, then proliferation, necrosis and cellular adhesion can balance to maintain spheroids and prevent further tumour fragmentation. Note that for sufficiently low levels of tumour aggressiveness (e.g., G = 0.10), tumour instability decreases until the steady-state configuration is a tumour spheroid, as predicted in  for non-necrotic tumours.
The finding that tumour morphology in the nutrient-poor regime may primarily depend upon the tumour microenvironment and not upon the tumour genetic/phenotypic characteristics has important implications for cancer treatment. In anti-angiogenic therapy, drugs are supplied to inhibit angiogenesis and the vascularization of the growing tumour and the host tissue in the microenvironment. Thus, anti-angiogenic treatment may result in hypoxia and heterogeneous nutrient delivery, effectively creating a nutrient-poor environment for the tumour. The resulting nutrient-poor microenvironment may induce tumour fragmentation, recurrence and metastasis. This is in fact observed in experiments, e.g. [57, 327]. This result is consistent with the findings of  who also found that anti-invasive therapy (increased adhesion, decreased mobility) may be used successfully as adjuvant to anti-angiogenic therapy.
The approach described in the previous section can be used to simulate tumour growth in a complex tissue. In figure 8, tumour growth in a heterogeneous domain that mimics brain tissue is considered (Macklin and Lowengrub ). In the white region (right side of the domain), μ = 0.0001, D = 0.0001 and B = 0 (the pre-existing blood vessel density), which models a rigid material such as the skull. In the black regions, μ = 10, D = 1 and B = 0, which models the cerebrospinal fluid. The light and dark grey regions model white and grey brain matter with regions μ = 1.5, D = 1 and B = 1 in the white and μ = 0.5, D = 1 and B = 1 in the dark grey. The tumour is denoted by a white thin boundary in the middle right of the frame. The proliferating, quiescent and necrotic regions in the tumour are coloured white, grey and black, respectively.
The simulations are from time t = 0 to t = 60 (approximately 45–90 days). The solution is plotted every 10 time units. The tumour grows rapidly until the nutrient level drops below nP = 0.30, at which time a large portion of the tumour becomes hypoxic and quiescent. The tumour continues to grow at a slower rate until the interior of the tumour becomes necrotic (see t = 10.0). This causes non-uniform volume loss within the tumour and contributes to morphological instability. Note that because the biomechanical responsiveness is continuous across the tumour boundary and the microenvironment has a moderate nutrient gradient, this simulation corresponds to the border between the invasive, fingering growth regime and the fragmenting growth regime that was described earlier. Additional effects can be seen that were not observed before, however.
As the tumour grows out of the biomechanically permissive tissue (light grey; μ = 1.5) and into the biomechanically resistant tissue (dark grey; μ = 0.5), its rate of invasion into the tissue slows (see T = 20.0). This results in preferential growth into the permissive (light grey) material, a trend which can be clearly seen from t = 30.0 onwards. When the tumour grows through the resistant tissue (dark grey) and reaches the fluid (black) (t = 40.0), the tumour experiences a sudden drop in biomechanical resistance to growth. As a result, the tumour grows rapidly and preferentially in the 1/2 mm fluid structures that separate the tissue (t = 50.0–60.0). Such growth patterns are not observed when simulating homogeneous tissues. Other observed differences are due to the treatment of quiescent (hypoxic) tumour cells. Regions that had previously been classified as necrotic (in [347–350]) are now treated as quiescent. As a result, tumour volume loss is reduced, and in particular, this may result in tumours with large hypoxic regions and little or no viable rim. Had these regions been treated as necrotic, the invasive fingers would have been thinner and the tumour may have fragmented. Therefore, the separate treatment of the hypoxic regions can have a significant impact on the details of invasive tumour morphologies.
The effects of extracellular matrix on tumour invasion were modelled by Anderson et al  by introducing a non-diffusing concentration field of bound matrix macromolecules E. The tumour cells degrade and remodel the ECM and move in response to gradients of E (haptotaxis). Along the same lines, Habib et al  employed a model describing chemotactic (motion in up gradients of soluble chemokines) and haptotactic cell behaviour, also without considering cell adhesion, and simulated tumour cell motility guided by the principle of least resistance, most permission and highest attraction. Earlier work by McElwain and Pettet  showed that in chemotaxis of cells in symmetric multicell spheroids may act against cell-motion induced by pressure gradients. Castro et al  developed a mathematical model of chemotactically directed tumour growth showing that heterotype chemotaxis provides an instability mechanism leading to the onset of tumour invasion, while homotype chemotaxis enhances the mean speed of the tumour surface, and tumour cell proliferation alone cannot generate the invasive branching observed experimentally. Khain et al  used a reaction–diffusion model to simulate chemotaxis and random cell motion of glioma. In agreement with recent experiments (Stein et al , unpublished results), they find that an outer invasive zone consisting of migrating cells grows faster than a higher-density inner region characterized by more proliferative cells. Their model predicts, however, that this is a transient state and the growth velocities of each region tend to the same value at long times. When the ratio of diffusion and cell-diffusion coefficients exceeds a critical value, Khain and Sander  find that symmetric fronts become unstable leading to instabilities and fingering as observed in experiments (e.g. [149, 152]). Chaplain and Lolas  developed a continuum model of ECM degradation by a particular matrix degrading enzyme (urokinase plasminogen) and studied the resulting tumour progression. Marchant et al  developed a model for tumour invasion and investigated the stability of travelling wave tumour–host fronts. In particular, the model predicts a biphasic dependence of the travelling wave speed with the density of the surrounding host tissue.
Propagating tumour–host wave fronts have also been investigated in the context of linear diffusional models of glioma tumour cells (e.g. Cruywagen et al , Woodward et al , Burgess et al , Swanson et al [508–512] and Jbabdi et al ). In particular, the effect of anisotropy on the diffusional propagation of glioma cell fronts has been investigated by Jbabdi et al (2005). Anisotropy may arise, for example, due to the presence of fibre tracks in the brain along which cells may preferentially migrate. Using diffusion tensor imaging (DTI) data of the brain, Jbabdi et al developed an anisotropic model for the cell-diffusion tensor. They find that the best fit between simulation results and clinical data occurs when the anisotropy of the cell-diffusion tensor is larger than the water diffusion tensor (which is directly imaged by DTI. Very recently, Swanson et al (2008) used patient data to parametrize the diffusion model in terms of the ratio D/ρ, where D is the net rate of dispersal and ρ is the net rate of proliferation. Swanson et al  then used images of untreated glioblastoma to provide initial conditions for the model and simulations of the subsequent growth were performed with and without resection (removal of the bulk of the tumour). A comparison between the actual and simulated tumours shows that the model is capable of accurately predicting patient survival (in the simulation, patient mortality is assumed to occur when the tumour radius is approximately 3 cm).
Using Lotka–Volterra type reaction–diffusion equations  that describe interactions between the tumour and host cells, Gatenby et al  assumed that the tumour–host interface is a travelling wave front and used an inverse problem approach to infer constraints on key biological quantities which appear as parameters in the model equations. In recent work, Zhu et al  proved the existence of travelling wave solutions to a system of equations originally proposed by Sherratt  and used later by Chaplain and Sherratt . In this approach, nonlinear diffusion is introduced to model contact inhibition between tumour and host tissues.
Boushaba et al  developed a two-compartment mathematical model to describe the suppression of secondary tumours nearby a large primary tumour. In their model, each tumour is assumed to stay in its own compartment and the compartments couple only through flux boundary conditions for the concentrations of growth and inhibitory factors expressed by the tumours, which diffuse through the extracellular matrix. Thus, reaction–diffusion equations describing the chemical species and tumour cell growth are posed in each compartment and flux boundary conditions are posed for the diffusing species. When a small secondary tumour is sufficiently close to the primary tumour, the inhibitory factor is found to suppress its growth. Removal of the primary tumour, and thus the source of the inhibitory factor, may result in the growth of the secondary tumour. If the secondary tumour is far from the primary tumour, the limited diffusion of the inhibitor prevents growth suppression and the secondary tumour is found to grow independently from the primary tumour.
Armstrong et al  recently constructed a continuous mathematical model of cell–cell adhesion by using non-local (integral) terms in a system of partial differential equations where cells use a so-called ‘sensing radius’ to detect their environment. Tindall and Please  analysed a mathematical model of avascular tumour spheroid growth that accounted for cell-cycle dynamics and chemotaxis-driven cell movement. Gerisch and Chaplain  applied a version of the model originally developed by Armstrong et al to formulate a continuum model of tumour cell invasion. This approach was further investigated by Szymanska et al . In , Macklin and Lowengrub developed a model to account for heterogeneous response to pressure and ECM adhesion gradients through non-constant cell mobility that depends on E and by introducing a haptotaxis velocity proportional to the gradient of E. For example, in equation (5) the cell velocity may be extended to read
Chemotaxis up a gradient of a soluble chemokine can be incorporated analogously (see also section 2.5.2). It is found that the morphology of the tumour is affected by heterogeneity in the ECM concentration, e.g. widespread variation in ECM density can result in the development of complex tumour morphologies. This is consistent with earlier findings by Anderson et al (2000). The ECM distribution also affects angiogenesis and vascular tumour growth (see section 2.4). Recently, Ambrosi  analysed cellular traction using an inverse problem approach. This was later used by Ambrosi et al  who applied a system of coupled elliptic partial differential equations to calculate the force field per unit surface generated by tumour cells on a polyacrylamide substrate. In the latter paper, the adjoint method is used. The shear stress thus obtained by numerical integration provided quantitative insight of the traction field and could enable study of the spatial pattern of force per unit surface generated in cell motion. Chauviere et al [117–119] developed models for cell migration due to chemotaxis and haptotaxis in anisotropic and heterogeneous tissues. In addition, using a growing and resting population model of cell migration was introduced to model the growth of glioma.
In recent work, Szymanska et al  developed a reaction–diffusion model to investigate the effect of heat shock proteins (HSPs) on cancer cell migration and solid tumour invasion. In vitro experiments were also performed. Cells that are placed in stressful conditions, e.g. elevated temperature and oxidative stress, tend to upregulate HSPs. HSPs perform multiple functions related to folding proteins into an appropriate shape and their upregulation in tumours is observed (e.g. ). HSPs are believed to act as a chaperone for many proteins linked to cancer progression including matrix-degrading enzymes (MDEs), and proteins involved in cytoskeleton remodelling. In experiments, inhibition of certain HSPs is found to reduce the rate of cancer invasion . Szymanska et al find that different mechanisms of HSP inhibition lead to different patterns of tumour invasion. Comparing these results with experiments suggests that HSP inhibition may decrease the cell motility without affecting the MDE production rate.
The effects of ECM and cell motility have also been incorporated in discrete tumour models (see section 3).
In hypoxic environments, tumour cells rely on anaerobic metabolism of glucose to lactic acid, even though this is much less efficient than aerobic metabolism. See Gatenby and Gawlinski [225, 227], Gatenby and Gillies  and Gillies et al . This can occur even under normoxic conditions. Although this is inefficient, it gives tumour cells a proliferative advantage over the host because tumour cells are typically able to survive in acidic environments that are toxic to the host cells. The first mathematical models used to study the effect of the glycolytic phenotype on tumour invasion potential were developed by Gatenby and Gawlinski  using continuum partial differential equations. In particular, a reaction–diffusion equation is introduced for the excess H+ ion concentration L, which assumes that excess ions are produced at a rate proportional to the tumour cell density and may diffuse and be reabsorbed:
where r is the acid production rate, NV is the tumour cell density, d is the reabsorption rate and DL is the H+ ion diffusion coefficient. Assuming as above that the tumour cell density is a constant, then NV is proportional to the characteristic function of the domain ΩV. The presence of acid in the microenvironment is assumed to cause death of both host and tumour cells, although at different rates. In the framework of the model described above, this can be introduced by appropriately modifying the tumour apoptosis parameter A to be dependent on L, and, analogously, introducing an apoptosis parameter for the host tissue. Relating L to the pH, Gatenby and Gawlinski  found that, consistent with experiments, the model predicts that gradients of acid develop in the microenvironment, which can trigger a transition to invasion. Further, an acellular gap may appear between the tumour and host tissues, as observed experimentally by Gatenby et al .
Smallbone et al  extended this model to incorporate a simple model of tumour vasculature. The vasculature is assumed to be a sink for H+ ions, which is found to prevent self-acidosis of tumours, thus suggesting that a possible treatment strategy to halt tumour growth is to introduce systemic acidosis. Gatenby et al  further extended the model to incorporate an extracellular matrix, which is degraded at the tumour–host interface by proteolytic enzymes released by acid-resistant tumour cells. This enables tumour cells to invade the damaged host tissue, consistent with experimental observations. In [487, 489] Smallbone et al incorporated the effect of tumour cell quiescence in the model. Since quiescent cells are less metabolically active, they are assumed to produce less acid than proliferating tumour cells. This may occur cyclically as tumour cells alternate from quiescent to proliferative states. The model is found to provide a description that better matches experimentally observed features of the tumour host interface, such as the size of the acellular gap between the tumour and host. There has also been modelling work to describe acidosis using discrete models (see section 3).
Heterogeneity in the tumour microenvironment may select for tumour cells that have a competitive advantage even if the underlying mutation rate is high, e.g. Nagy . Gatenby and Frieden  demonstrated this using information theory and showed that this is required in order to prevent the accumulation of genetic malfunctions. Gatenby and Vincent  developed a model using population biology and game theory to describe competing clonal populations that interact with one another and the microenvironment. The model further confirms that tumour invasion requires a certain degree of genetic stability. In particular, it is found that malignant tumours are more homogeneous genetically than their pre-malignant counterparts. Gillies and Gatenby  went on to propose that the evolutionary dynamics governing the microenvironment selection forces could be analysed using an inverse problem approach. In particular, the evolution dynamics can be inferred from the typical phenotypic traits of the tumour. This analysis suggests that adaptation to hypoxia and acidosis must be a major component of the carcinogenic sequence. It should also be noted that therapeutic intervention can also exert selection pressure on a multi-clonal tumour cell population.
Experiments demonstrate that mechanical stress may significantly influence the growth of solid tumours (e.g. ). The mathematical modelling of residual stress development in growing tumours was recently reviewed by Araujo and McElwain , Roose et al  and Tracqui . Briefly, the study of Shannon and Rubinsky  showed that in a linear-elastic description of a growing tissue with spherical geometry, such as a tumour spheroid in vitro, residual stresses are induced by any spatial variation in the growth process. Nonlinear elastic models were employed by Chaplain and Sleeman  to classify solid tumours and by Chen et al  to investigate the influence of the growth-induced stress in the medium surrounding a tumour spheroid on the growth of the tumour. The growth rate and equilibrium size of the tumour were found to decrease as the stiffness of the surrounding medium increased consistent with experiments (e.g. ).
Jones et al  considered the continuous nature of the growth process instead of a fixed growth strain distribution, extending the work of Shannon and Rubinsky. The dimensionless model equations from Jones et al (2000) are
where the first equation describes nutrient diffusion/uptake, the second describes the rate of volume growth, the third equation describes mechanical equilibrium where σ(x, t) is the stress tensor. The last equation is a differentiated version of the stress–strain relation with and I the identity tensor. The Jaumann derivative has been used to maintain frame-indifference (D/Dt = t + v · represents the material derivative). As before, n(x, t) is nutrient concentration and v(x, t) is cell velocity; the tumour–host interface Σ is advected with the velocity v. At the tumour/host interface traction-free boundary conditions are applied (although one could also apply a normal stress on the boundary to mimic the effect of a gel containing the tumour).
One of the difficulties with the model described above is that the stress does not evolve to a steady-state distribution even when the tumour reaches a steady shape. Incorporating mechanisms of stress relaxation can allow steady-state stress distributions to be achieved, however. This can be done by using poroelasticity models (see Roose et al ), by incorporating anisotropy (see Araujo and McElwain [35, 38]) or viscoelasticity (e.g. MacArthur and Please ). The effect of solid stress on the vasculature was investigated by Sarntinoranont et al  using the poroelastic model and by Araujo and McElwain  using the anisotropic model.
A significant difficulty in using the approach described above is that it relies on frame-invariant differentiation of the stress–strain relation. However, as Ambrosi and Preziosi  point out, the strain is not frame invariant  and thus the Jaumann derivative is inappropriate. Instead, the use of accretive forces (e.g., [13, 19–22, 24, 221, 279, 280, 320, 321, 340]) is required to derive the equations for growth in a consistent manner from a reference configuration. Interestingly, as pointed out by Ambrosi and Preziosi , the above models (with the convection terms dropped) may be obtained from a small deformation limit of the constitutive laws derived in  that take into account accretive forces.
In a single-phase tumour model, Ambrosi and Mollica [21, 22] utilized a nonlinear approach with accretive forces to analyse the role of mechanical stress on the growth of a tumour spheroid. A mechanical description was used where the volumetric growth and mechanical responses were divided into two separate contributions. In particular, the deformation tensor F is decomposed multiplicatively to account for the contributions of pure growth and elastic deformation:
where Fr describes the deformation due to cell reorganization and stress relaxation after growth, G describes the deformation due to growth. A Blatz–Ko constitutive equation  was used to relate the stress and the elastic part of the deformation Fr. The growth tensor was taken to be G = gI. Ambrosi and Mollica  derived the following equations for a tumour spheroid
where ρ is the tumour density, ρ0 is the tumour density at time t = 0, Jr = det(Fr), J = det(F) and P is the first Piola–Kirkhoff stress tensor obtained from the Blatz–Ko constitutive law. The derivatives in the above are with respect to the reference configuration, and the overdot denotes the time derivative. The cell velocity is v = , where u is the displacement vector. Normal stress boundary conditions may be applied at the tumour boundary, or the tumour may be embedded in an elastic matrix material.
Ambrosi and Mollica  applied this system to describe homogeneous growth inside a rigid cylinder, modelling ductal carcinoma, and to the growth of a tumour spheroid with non-homogeneous diffusion of nutrients, which generates residual stresses because the nonuniform distribution of nutrients leads to inhomogeneous growth. In later work , these authors provided a qualitative analysis of the stress-modulated growth of a continuum body as predicted by equations that satisfy an a priori stated dissipation principle. Accretive forces  were re-interpreted as a homeostatic value of the Eshelby stress, coinciding with the classical biomechanical concept of homeostatic stress in the case of infinitesimal strain.
The previously reviewed models were all in one-dimension, i.e. a spherically symmetric geometry. A general stability analysis of nonlinear elastic models to asymmetric perturbations was performed by Ben Amar and Goriely  in higher dimensions. The important role of residual stress was established by showing that a spherical shell without any external loading could become spontaneously unstable under large anisotropic growth.
In recent work, Lloyd et al  developed a finite-element based method to simulate the elastic response of tissue during tumour growth in 2D. In particular, deformations are induced by a prescribed strain determined from cell proliferation. However, the accumulation of stress is neglected and it is assumed that residual stress dissipates completely, which would be consistent with an elasto-plastic growth law.
A number of other models of stress effects in tumours have been developed in the context of multiphase mixture models; see section 2.5.
To transition from the avascular to the vascular phase of growth, a tumour must induce new blood vessels to sprout from the existing vascular network and grow towards the tumour, eventually penetrating it. This process, known as tumour-induced angiogenesis, is a critical milestone in the development of invasive and malignant cancer . The process is thought to start when a small avascular tumour exceeds a critical size greater than can be sustained by the normal tissue vasculature . Accordingly, tumour cells become hypoxic and secrete diffusible chemical factors, such as the vascular endothelial cell growth factor (VEGF). The factors diffuse into the host microenvironment and bind to specific membrane receptors on the (vascular) endothelial cells (ECs) that line existing blood vessels. This activates migration of ECs, which respond by degrading the basement membrane surrounding the existing vessel to form sprouts. The ECs then proliferate and migrate towards the tumour. Migration is mediated by the chemotactic response to VEGF and other pro-angiogenic factors, by proteolytic enzymes that degrade the ECM providing space for the cells to move, and a haptotactic response to variable cell–matrix adhesion. As the ECs proliferate through the ECM, they form tubular structures which fuse (anastomose) to form loops. Eventually blood flows through the neovascular network providing the tumour and host microenvironment with an increased, although inefficiently delivered, supply of cell substrates. This process is illustrated in figure 9.
Tumour growth and angiogenesis are coupled in that hypoxic tumour cells release a net balance of pro-angiogenic factors that attract the ECs and incite the neovascular network to approach the tumour. This creates additional sources of cell substrates in the microenvironment. The tumour responds by upregulating cell proliferation in regions where nutrient is increased. The additional oxygen and nutrient affect the hypoxic tumour regions, which in turn downregulates the net release of pro-angiogenic factors and hence affects the formation of the neovascular network.
In addition to heterogeneous blood flow, a neovascular network contends with pressure variation introduced by increased tumour cell proliferation and migration. Normally, the network responds by remodelling itself. Compared with vessel networks formed during normal biological processes during development and wound healing, however, tumour-induced neovascular networks may become leaky and inefficient, producing immature and tortuous vessels , which leads to increased flow resistance and ultimately to heterogeneous supply of cell substrates in the tumour microenvironment .
Mathematical models of tumour-induced angiogenesis date to the work of Balding and McElwain . Both continuum, fully discrete, composite continuum–discrete and hybrid continuum–discrete mathematical models have been developed (see recent reviews [15, 113, 299, 332, 333, 362, 406, 410, 427, 429]). In most models, the coupling between tumour growth and angiogenesis is simplified in that one of the two processes is static while the other is dynamic. There are few models in which both processes are coupled dynamically.
Models that focus on the angiogenic response have taken two approaches. One approach focuses on blood vessel densities rather than vessel morphology, as in continuum partial differential equations. See for example Byrne and Chaplain , Orme and Chaplain , De Angelis and Preziosi , Sansone et al , Levine et al [334, 335, 337], Hogea et al , Peterson et al , Stamper et al , Jain et al  and Addison-Smith et al . In these approaches, continuum conservation laws are introduced to describe the dynamics of the vessel densities and angiogenic factors. These models do not provide morphological or blood flow information about the vasculature. In the context of vasculogenesis in vitro, biomechanical and biochemical models have been developed that account for cell–ECM interactions and chemotactic response and are capable of describing the morphology of the vasculature. In these models, the network morphology emerges, roughly speaking, from a homogeneous distribution through a type of phase transition. See, for example, [17, 125, 217, 276, 329, 358, 386, 392, 475, 524] for continuum models and [373, 374, 375] for discrete models.
The other approach represents vessels as line segments, continuous curves, interconnected lattice patterns or collections of individual endothelial cells. Mechanisms have been modelled such as vessel sprout branching and anastomosis, as well as vascular endothelial cell activation, proliferation and migration via chemotaxis up gradients of tumour angiogenic factors (e.g. VEGF), haptotaxis up gradients of ECM-bound chemokines (e.g. fibronectin) and proteolysis of the ECM. See, for example, Stokes and Lauffenberger , Anderson and Chaplain , Tong and Yuan , Plank and Sleeman [417, 418], Sun et al [502, 503], Kevrekidis et al , Bauer et al , Milde et al  and Capasso and Morale . Blood flow and network remodelling have also been simulated (e.g. Pries et al , McDougall et al , Stephanou et al [494, 495], McDougall et al , Wu et al , Zhao et al , Sun and Munn  and Pries and Secomb ). In addition, models of tumour growth in static network topologies have been simulated (e.g. Alarcón et al , Betteridge et al ) and multiscale models of fluid transport in tumour have also been developed (Chapman et al ).
The first model in which tumour growth is fully coupled with tumour-induced angiogenesis for an arbitrary network topology was presented in Zheng et al . This was modelled in 2D using a version of the continuum–discrete model of Anderson and Chaplain  coupled to the non-symmetric tumour growth model of Cristini et al . This approach is discussed in further detail below. Later, in the context of discrete cell-based systems, models of angiogenesis and vascular tumour growth were also implemented, see, for example, Gevertz and Torquato , Bartha and Rieger , Bauer et al  and Scislo and Dzwinel . The effects of blood flow through the neovascular network on tumour growth were also recently considered in Alarcón et al , Bartha and Rieger (2006), Lee et al , Welter  and Owen et al  using cellular automaton models for tumour growth combined with dynamic network models for the vasculature. The effects of an arterio-venous network were considered in Welter et al . Very recently, Macklin et al , extended the model of Zheng et al (2005) and incorporated a version of the dynamic model of tumour-induced angiogenesis developed by McDougall et al  to explicitly analyse the effects of blood flow and vessel remodelling. In other recent work, Lloyd et al  simulated the vascular growth of a 2D tumour by coupling models for angiogenesis, flow through the developing neovascular network and network remodelling with an elastic tumour growth model they developed earlier (Lloyd et al ). In 3D, vascular tumour growth has been studied using a mixture model and lattice-free description of tumour-induced angiogenesis (Frieboes et al  and Bearer et al ).
A basic angiogenesis model can be constructed  based on that of Anderson and Chaplain  with some additions taken from Paweletz and Knieriem , Paku , Chaplain and Stuart ( ) and McDougall et al . The model of Anderson and Chaplain  uses a composite continuum/discrete approach with the ability to follow the motion of individual endothelial cells at the capillary tips and control important processes such as migration, proliferation, branching and anastomosis using a discrete random walk algorithm. Cell substrates such as tumour angiogenic factors and ECM are described using continuum fields. Their objective was to replicate angiogenesis as observed in the experimental ‘rabbit eye model’ , where a tumour is implanted in the cornea, thus inducing angiogenesis that can be readily observed (since the cornea is normally avascular). This is artificial because most tissue in the body is normally well vascularized.
Accordingly, the following field variables were introduced: 
Once a tumour cell senses that the nutrient level has dropped below a critical threshold, the cell releases diffusible tumour angiogenic factors (TAFs). TAF may be described through a reaction–diffusion equation with a point or line boundary condition at the necrotic/viable tumour cell interface. The TAF molecules are much smaller than cells and diffuse quickly through the extracellular spaces. As a result, a quasi-steady reaction–diffusion equation can be assumed for the non-dimensional concentration T(x, t) of TAF 
where is the diffusion coefficient, 1sprout tips is the characteristic function of the sprout tips, denote the non-dimensional production, natural decay and binding rates of TAF, and is the production rate of TAF by the tissue (see equation (42)). In the far-field at the boundary of the computational domain, zero Neumann boundary conditions T/n = 0 may be taken.
A primary component of the extracellular matrix is fibronectin, a long, non-diffusible binding molecule. ECs produce, degrade and attach to these molecules during their migration towards the tumour. The concentration of fibronectin f(x, t) satisfies 
where ηp is rate of production of fibronectin by ECs and ηm is the rate of degradation of fibronectin by MDEs . The MDE concentration, m, satisfies
where the terms on the right-hand side represent diffusion (MDEs diffuse slowly and so the quasi-steady approximation is not appropriate), production of MDE by viable tumour cells and sprout tips and natural decay of MDE, respectively.
While ECs are comparable in size to host cells and tumour cells, it may be assumed  that there are not enough ECs to modify the cell velocity u. The ratio of endothelial to tissue cells is of the order 1/50 or 1/100 .
At the continuum level, the density of endothelial cells obeys a reaction–diffusion–convection equation. The problem is convection-dominated, with primary source of convection driven by chemotaxis and haptotaxis of endothelial cells in response to TAF and ECM (fibronectin) concentrations, respectively. In addition, the ECs may also be affected by the cell velocity. At the continuum level, the density e(x, t) of ECs, which is related to the probability of finding the tip of a capillary at that location and time, obeys a convection–reaction–diffusion equation in ΩH and ΩV :
where De is the EC diffusion coefficient, χc, χf are chemotaxis and haptotaxis coefficients, respectively, χu and α are dimensionless coefficients and is the Heaviside function. The parameter χu measures the degree to which the capillaries are influenced by the cell velocity. Further, c* is a concentration of TAF above which proliferation occurs, ρD, ρp, and ρN are rates of natural degradation, production and necrosis of ECs, respectively.
In the original work of Anderson and Chaplain , the model for motion of capillary sprout tips comprised continuum and discrete components. Equations (36), (37) and (39) constitute the continuum component. The discrete component was derived from equation (39) under the assumption that growth of the capillary is determined by the biased random migration (random walk) of a single endothelial cell at the sprout-tip. In particular, it was assumed that there is a trail of ECs that follow the sprout-tip. In later work, McDougall et al  omit the continuum EC equation (39) and simply use the discrete approach. In the discrete algorithm, probabilities are generated from a finite-difference approximation that describe the tendency of the tip endothelial cell to migrate and proliferate on a Cartesian lattice. In addition, capillary branching and anastomosis are incorporated. The entire capillary may be convected by the external cell velocity using the kinematic condition 
where x is the position on the capillary and μC is the capillary mobility . As the vessel network becomes more established in time, the capillaries mature into larger vessels and become more rigid, which can be accounted by decreasing μC.
Alternative and related approaches have also been developed. Plank and Sleeman [417, 418] and Plank et al  considered lattice-based and lattice-free reinforced random walk models of angiogenesis. Sun et al [502, 503]) coupled deterministic models of ECs and vessels at the cell-scale with partial differential equation models of chemical factors at the tissue-scale. Kevrekidis et al  and Whitaker et al  incorporated the effects of inhibitors, following earlier work by Levine , in continuum–discrete models of angiogenesis. Bauer et al  modelled angiogenesis using the Graner–Glazier–Hogeweg (GGH) framework of cell-level modelling (see section 3). Milde et al  incorporated the effects of both matrix-bound and soluble growth factors in a three-dimensional model of angiogenesis. More recently, Capasso and Morale  developed a multiscale framework for modelling angiogenesis that couples discrete models with stochastic effects at the cell-level and coarse-grained partial differential equation models at the tissue scale.
The flow in the neovascular network and its effects on network remodelling have been explicitly simulated using network models by Pries et al , Pries et al [432, 433]. Following this work, McDougall et al  extended the basic angiogenesis model above by incorporating blood flow by treating the network as a series of straight, rigid cylindrical capillaries that join adjacent nodes. Blood flow was modelled through the cylindrical vascular network by representing the elemental flow-rate in each segment with Poiseuille's Law, which describes flow rate as a function of capillary lumen, fluid viscosity, capillary length and pressure drop. Stephanou et al  extended this model to three dimensions and found that the highly interconnected nature of irregular vasculature produced by tumour-induced angiogenesis could cause low rates of blood flow to the tumour with the potential for blood-borne drugs to bypass the entire mass depending on the tumour shape. They also examined the effect of vessel pruning (e.g. as may occur during anti-angiogenic therapy) on flow through the vascular network. In order to investigate how adaptive remodelling affects oxygen and drug supply to tumours, these authors later included vascular adaptation effects  (due to shear and circumferential stresses generated by flowing blood [433, 434]). This model was further updated by McDougall et al  to incorporate dynamic vessel radius adaptation, thus coupling vessel growth with blood flow, in contrast to earlier flow models where effects of blood flow were evaluated after generating a hollow network (e.g. as in Stephanou et al  and Alarcón et al ).
The nonlinear coupling of vessel growth and adaptation with flow has a nontrivial effect on the vascular morphology and development. In order to couple angiogenesis and tumour growth, the nutrient transfer from the vascular system to the host and tumour tissues needs to be modelled. In Zheng et al  and Macklin et al , angiogenesis models are coupled with the tumour growth model described in section 2.3. The nutrient equation (20) is modified to account for nutrient released by the neovascular network. Accordingly, following Macklin et al, an additional source term is added to the right-hand side of equation (20)
where is the transfer rate from the neovascular vessels. The function Bneo (x, t) = 1neo is the characteristic or indicator function of the neovasculature (i.e. equal to 1 at the locations of the new vessels). Further, P is the oncotic (solid/mechanical/hydrostatic) pressure, Pvessel and h are the dimensional pressure and the haematocrit in the neovascular network, respectively. The constants D and min reflect the normal value of haematocrit in the blood (generally about 0.45) and the minimum haematocrit needed to extravasate oxygen, respectively. The haematocrit may be 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 (Pvessel, P) such that large oncotic pressures P relative to vessel pressures Pvessel may prevent extravasation and transfer of oxygen from the vessels into the tissue.
The oxygen source term in equation (41) is designed such that for a sufficiently large transfer rate , the oxygen concentration n ≈ 1 at the spatial locations of the neo-vessels. Note that oxygen flux conditions across the neovasculature could be imposed instead, see, for example, Alarcón et al .
Another form of coupling between the growing tumour and developing vascular network occurs through the angiogenic factors released in the tumour microenvironment. In Zheng et al  and Macklin et al , hypoxic/quiescent tumour cells are assumed to secrete tumour angiogenic factors (TAFs), which diffuse into the surrounding tissue and attract endothelial cells (ECs). ECs respond by binding with the TAF, proliferating and chemotaxing up the TAF gradient. Following Macklin et al  the source term in equation (36) describes the source of TAF due to secretion by hypoxic cells in the perinecrotic region:
where 1ΩQ is the characteristic function of the region of quiescent (hypoxic) cells and denotes the non-dimensional production rate of TAF.
The effect of solid/mechanical pressure-induced vascular response on tumour-induced angiogenesis and vascular growth is shown in figures 10 and and1111 from Macklin et al . With the values of the parameters used here, a solid pressure-induced vascular constriction occurs when the pressure P ≈ 0.8. Angiogenesis is initiated from an avascular tumour configuration at t = 45 days (not shown), when 10 sprout tips are released from the parent vessel. At early times, 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 (55 days of total growth time). Blood flow in the neovasculature starts near the parent capillary and eventually the flow reaches the tumour.
The solid pressure prevents delivery of oxygen internally to the tumour, and thus the delivery of oxygen is heterogeneous and significant 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 neovasculature, the solid pressure generated by tumour cell proliferation constricts the neovasculature in the direction of growth (where pressure is highest) and also correspondingly inhibits the transfer of oxygen from those vessels. This makes the tumour grow more slowly.
The neovasculature in other areas of the host microenvironment then provides 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. 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 complex tumour morphologies. This result is consistent with the theory and predictions made earlier that substrate inhomogeneities in the tumour microenvironment tend to cause morphological instabilities in growing tumours. See, for example, Byrne and Chaplain [89, 90], Cristini et al [126, 129, 339], Anderson et al [26, 32, 235] and Macklin and Lowengrub [348–351].
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 the nearby neovasculature in the host microenvironment. The neovascular response is 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 figure 11, the dimensional neovasculature radii (in m) and intravascular pressures (in Pa) are shown together with the non-dimensional ECM and TAF concentrations. Blood flow causes a dilation of the vessels and an overall decrease in pressure as branching, anastomosis and increased blood flow occur throughout the neovascular 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 its interior. The new vessels are still able to migrate through the region of lower ECM even though this acts against haptotaxis. Because the tumour grows slowly, only the tips of the invasive clusters outrun the degraded ECM. As can be seen in figure 11, 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 strong nonlinear coupling between the tumour-induced angiogenesis and the progression of the tumour. The pressure-induced vascular response of constricting the radii of the neovasculature and inhibiting blood-tissue oxygen transfer not only affects the tumour growth dramatically, but also significantly affects the growth of the neovascular network.
The version of this model developed by Zheng et al (2005) has also been used by Cristini et al  to predict changes in tumour morphology in response to perturbations in the model parameters that govern cell–cell adhesion and the density of the microvasculature in the host tissue, leading to the creation of a ‘morphology diagram’  (figure 12). If cell–cell adhesion is sufficiently strong (case A), then the tumour tends to maintain a compact morphology, even following the onset of angiogenesis. In contrast, when cell–cell adhesion is low (case B), the tumour tends to break into fragments  that invade the surrounding tissue due to cell substrate gradients [93, 484, 575]. When adhesion is kept low but the host vascular density is increased (case C), the substrate levels become more homogeneous (were ‘normalized’), leading to more compact morphologies and reduced invasion. These parameter ranges can be achieved through therapeutic intervention. In anti-invasion therapy, drugs are introduced to increase cell–cell adhesion; in anti-angiogenic therapy, drugs are used to target and interrupt the tumour neovasculature; and in vascular-normalization therapy, inefficient blood vessels are ‘pruned’ to reduce or eliminate hypoxic gradients. These results may help explain some undesirable effects of current anti-angiogenic therapies, due to exacerbating hypoxic gradients . For example, the model predicts that anti-invasive therapy will lead from a transition from either case B or C to case A, while anti-angiogenic therapy causes a transition from case C to case B, and vascular-normalizing therapy reverses the transition.
Following the strategy of Zheng et al , Frieboes et al  simulated angiogenesis and vascular tumour growth in 3D using a multiphase tumour model  coupled with a lattice-free continuous–discrete model of angiogenesis originally developed by Plank and Sleeman . The angiogenesis model involves a random walk on the unit sphere (instead of on a Cartesian lattice). The lattice-free angiogenesis model is otherwise similar to the basic angiogenesis model described above.
The models reviewed so far treat the tumour mass as a single-component material that locally expands and contracts in correspondence to variable rates of cell adhesion, mitosis, necrosis and apoptosis. In fact, tumours consist of multiple phases including a variety of different cell genotypes and phenotypes as well as ECM and water. Being able to describe the dynamics of multiple cell species is crucial since the tumour microenvironment (e.g. hypoxia) and impaired cell genetic mechanisms can lead to multiple cell genotypes and phenotypes that select for cell survival under abnormal conditions , with profound consequences on the overall tumour growth, invasion and response to treatment. Furthermore, the microenvironment of invasive tumours may be characterized by non-sharp boundaries between tumour and host tissues, and between multiple species within the tumour [323, 327, 341]. A multiphase approach represents a more general and natural modelling framework for studying solid tumour growth, and is able to provide a more detailed account of the biophysical process of tumour growth than that in single-phase models. Thus, this is an important direction for future research.
In a multiphase approach, a solid tumour is described as a saturated medium, comprising of at least one solid phase (cells, ECM, etc) and one liquid phase (water), and can be generalized to incorporate any number of additional phases to describe multiple cell species. At a given location, the mass or volume fractions describe the relative amounts of different cell clones as well as necrotic and host tissue. The governing equations consist of mass and momentum balance equations for each phase, interphase mass and momentum exchange, and appropriate constitutive laws to close the model equations. Phenotypic and genetic heterogeneity (e.g. mutations) in the cell population can be represented. This is a critical improvement towards realistically simulating mutation-driven heterogeneity. The mixture model approach also eliminates the need to enforce complicated boundary conditions across the tumour/host (and other species/species interfaces) that would have to be satisfied if the interfaces were assumed sharp. Further, this methodology does not require explicit tracking of the interface as is required with a sharp interface.
Recently, multiphase mixture models have been developed to account for heterogeneities in cell-type and in the mechanical response of the cellular and liquid tumour phases. In early work, Please et al [420, 421] applied multiphase modelling to tumour growth by capturing both tumour cells and extracellular fluid as separate continuum phases. Ward and King [543, 544] and Breward et al  modelled avascular cancer growth through a two-phase description, comprised of tumour tissue and dead tissue (extracellular space) and incorporated a model of cell–cell adhesion. This model was later extended by Venkatasubramanian et al  and Bertuzzi et al  to account for ATP production from energy metabolism in the models of multicellular tumour spheroids.
Ambrosi and Preziosi  reviewed several mixture models and developed a mixture model treating the tumour as a deformable porous material. Byrne et al  and Byrne and Preziosi  extended the two-phase model of Ambrosi and Preziosi (2002). Jackson and Byrne  and Lubkin and Jackson  used incorporated the biomechanical effects of a capsule surrounding a tumour and analysed the effects of encapsulation on tumour growth. Jackson and Byrne (2002) also investigated the consequences of different hypotheses for capsule formation. Byrne and Preziosi (2003) investigated the effect of stress-dependent cell proliferation and external loads on spherical tumour growth and the equilibrium size attained. Franks et al [186, 187] used mixture models to study ductal carcinoma in situ of the breast. Breward et al  extended avascular modelling to describe vascular tumour growth, thus incorporating a third phase to describe the spatial and temporal distribution of blood vessels.
Roose et al  studied the stress generated by solid tumour growth using a poroelastic model. Araujo and McElwain [37, 38] proposed an alternative multiphase model of tumour growth in an effort to more accurately capture residual stresses using an energetically derived thermodynamically consistent mixture model that includes a solid phase representing the tumour cells and extracellular matrix and a liquid phase. Isotropic and anisotropic growth models were considered, highlighting the need to incorporate stress relaxation in order to predict a stable evolution of stresses over a period of growth and equilibration to a steady avascular state. As discussed earlier, however, these methods rely on an inappropriate use of frame-indifference and the Jaumann derivative and instead accretive forces should be used . Chaplain et al  used a mixture model to focus on biomechanical effects on tumour growth and demonstrated that if tumour cells underestimate the compression state of the tissue, this may confer a fitness advantage.
Bertuzzi et al [63, 64] developed limited multiphase models of tumour cords (growth of a tumour around a capillary) that incorporate cytotoxic agents and the flow of an interstitial fluid. Geometry-dependent kinematic relations limit the generality of these models. Astanin and Tosin  developed a fully multiphase mixture model of tumour cords growing around a capillary. Later, Preziosi and Tosin  developed a model of adhesive interactions between tumour cells and the ECM and performed two-dimensional simulations to investigate the effect of these interactions on the growth of tumour cords and the development of fibrosis. Astanin and Preziosi  used a two population model to investigate the upregulation of glycolysis (Warburg effect) in tumour cords. Ambrosi and Preziosi  used a mixture model to study the elasto-viscoplastic response of the tumour and microenvironment during growth. Preziosi et al  extended previous work and developed an elasto-visco-plastic model of cell growth. This was then used to compare with experimental results of cell aggregation.
Green et al  developed and analysed a multiphase mixture model to investigate the effect of interactions among cells and the ECM on the formation and structure of multicellular aggregates in vitro. Galle et al  used a multiphase mixture model to study the effect of contact inhibition on the growth of cell colonies and compared the results with an agent-based discrete model of cell growth. Further references on this type of modelling can be found in the recent reviews by Araujo and McElwain , Hatzikirou et al , Quaranta et al , Byrne et al , Graziano and Preziosi , Astanin and Preziosi , Roose et al , Preziosi and Tosin  and Tracqui . Due to the complexity of multiphase models, most analyses and numerical simulations are one-dimensional or radially symmetric.
Very recently, thermodynamically consistent mixture models for all phases of solid tumour growth in 3D have been developed. Wise et al  developed a multispecies mixture model and performed simulations of tumour growth in three dimensions. Frieboes et al  and Bearer et al  combined the mixture model developed in Wise et al with a model for tumour-induced angiogenesis to simulate vascular growth and compared results with clinical data. Cristini et al developed, analysed and simulated a mixture model to study tumour invasion and branching under hypoxic conditions. Using a general approach based on energy variation, the nonlinear effects of cell-to-cell adhesion and taxis inducing chemical and molecular species have been incorporated. This model enables a detailed description of tumour progression and the dependence of cell–cell and cell–matrix adhesion on cell phenotype and genotype as well as on the local microenvironmental conditions (e.g. oxygen levels). The system energy accounts for all the processes to be modelled. Adhesion is introduced through an interaction energy that leads to well-posed fourth-order equations. In this approach, sharp interfaces are replaced by narrow transition layers that arise due to differential adhesive forces among the cell species. The resulting diffuse-interface mixture equations are well posed, unlike some previous mixture models, and are a coupled system of equations including fourth-order nonlinear advection–reaction–diffusion equations of Cahn–Hilliard-type  for the cell-species volume fractions coupled with reaction–diffusion equations for the substrate components. A related non-local continuum model of adhesion was recently developed by Armstrong et al , and used later in the tumour context by Gerisch and Chaplain . We also note that a continuum model of cell adhesion and migration was recently developed by Kuusela and Alt  in the context of a slowly migrating cell on a substrate.
To account for angiogenesis, this mixture model has been coupled nonlinearly to a composite continuum–discrete, lattice-free model of tumour-induced angiogenesis originally developed by Plank and Sleeman [417, 418]. To solve the equations numerically, very efficient, adaptive finite-difference nonlinear multigrid methods have been developed [128, 553–555]. The model was employed to study 3D vascularized cancer growth (Frieboes et al ), including malignant brain tumours (Frieboes et al , see also Sanga et al ).
The primary dependent variables in a (N + 1)-species model can be specified as follows :
It may be assumed that there are no voids (i.e. the mixture is saturated) and thus . Further, for simplicity, it may be assumed that densities are constant and equal to ρ, i.e. independent of temperature, pressure and component type. We also suppose that the system is isothermal. Without loss of generality, i = 0 is identified as the water component. The volume fractions of the components are assumed to be continuous in a domain Ω, which contains both the tumour and host tissues.
The volume fractions obey the mass conservation (advection–reaction–diffusion) equations
where Ji are fluxes that account for the mechanical interactions and diffusional fluxes among the cell species. The source terms Si account for inter-component mass exchange as well as gains due to proliferation of cells and loss due to cell death and lysing.
The volume-averaged velocity of the mixture is then defined as . Summing equation (43), the mass of the mixture is conserved only if is constant. Without loss of generality, this constant is chosen to be zero:
These conditions are posed as consistency constraints for the fluxes and sources. In the absence of inertial and external forces, the balance of linear momentum is
where πi are interaction forces among the species. Assuming that the mixture stress satisfies · σ = 0, this gives the constraint .
Let ui be the internal energy of the ith component and let the volume-averaged internal energy of the mixture be . Then, in its simplest form, the balance of energy equation is given by (e.g. [37, 128])
where Di/Dt = (/t) + ui · is the advective derivative with respect to the velocity ui, ri is the heat added to each phase to keep the mixture isothermal (θ is the fixed temperature) and βi are energy interaction terms that satisfy
Constitutive relations for the fluxes Ji, the interaction forces πi and the interaction energies βi may be posed consistently with the second law of thermodynamics. Briefly, the idea is as follows. Let ηi denote the entropy of each component. Then, the volume-averaged entropy of the mixture is . Defining the temperature to be θ (which is assumed to be constant), the second law of thermodynamics can be posed in terms of the Clausius–Duhem inequality 
where is an entropy flux and , where ri are the volume-averaged rates of heat supply needed to keep each component isothermal. The entropy and internal energy ui may be used to define the Helmholtz free energy of each component ψi = ui – θηi and the free energy per unit volume Ψi = ρϕiψi. The next step is to re-write the energy equation (46) in terms of the Helmholtz free energy and to posit the dependence of the free energy upon independent variables of each phase .
For binary mixtures of one solid (i = 1) cellular component and water (i = 0), Araujo and McElwain , assume that ψi = ψi(F1, G1, u1 – u0). In this approach, F1 is the deformation gradient in the solid phase and G1 = F1. If one further assumes that the volume fractions of solid and liquid phases are constant, Ji = 0, and that the liquid phase is inviscid, then Araujo and McElwain show that the stresses become
where p is a pressure that is introduced to maintain ϕ0 + ϕ1 = 1 and I is the identity tensor. Assuming that displacements are small, and that the solid component is incompressible, Araujo and McElwain went on to derive the mixture system:
for conservation of mass. The growth equations are
where is the Jaumann derivative (e.g. ), Ω is an anisotropic growth tensor  and μ is a Llame constant. The stress–strain relation is E1 = (1/2μ)σ1−((trσ1 + 3ϕ1p)/6μ) I + gΩ, where g is the rate of volume growth of the solid which satisfies
Finally, the momentum equations are
where α is a drag coefficient.
An alternative approach was taken by Cristini et al . In their work, the solid fraction was treated as a viscous fluid and cell–cell and cell–matrix adhesive interactions were incorporated and thermodynamically and mechanically consistent equations were derived. In addition, taxis inducing chemical and molecular species were included. Accordingly, the Helmholtz free energy was taken to be
where c1, . . . , cL are the concentrations of chemical and molecular species, χil are taxis coefficients and εil measures the strength of component interactions (see Wise et al ). The dependence of Ψi upon the volume fractions and the gradients of the volume fractions arises naturally through expansion of a nonlocal interaction potential among the phases (e.g. see [98, 555]) that represents adhesive interactions among the species. The resulting system is
and . The flux in the liquid is . In the liquid (assumed to be inviscid), the momentum balance equation is the multicomponent Darcy's law:
where the αi are drag coefficients. For the solid components, the general viscous law is obtained (for i > 0):
where , where λi and vi are viscosities.
In the case of two component mixtures, a local approximation of this model may be achieved by taking u0 = −(ϕ1/ϕ0)u1, p to be a constant, setting the viscosity to zero and dropping equation (56) to get u1 = −Mϕiμi, see . Viscosity can be easily included and yields a nonlocal equation analogous to the mixture model derived in Byrne and Preziosi  in one-dimension; the inviscid model is similar to that considered earlier by DeAngelis and Preziosi . Interestingly, the local approximation can be shown to converge to a classical single-phase sharp-interface model of the type considered in previous sections as the component interactions εil = ε → 0, see Cristini et al . A discussion of the model without the above approximation may also be found in .
For multicomponent mixtures consisting of more than two components, one may take yet another approximation by supposing that the velocities of the solid components are given by ui = Miϕiμi and that the velocity of the liquid is . This is analogous to the closure laws described in DeAngelis and Preziosi (2000) and Ambrosi and Preziosi .
Another approach to modelling multicomponent tumours was recently developed by Frieboes et al and Wise et al [193, 555]. In this approach, the mass balance equations are the same as in equation (54) and the constitutive laws for the mechanical fluxes Ji and a generalized Darcy's law for the cell velocities ui are derived using an energy variation argument utilizing the mixture energy E given by
and , where i > 0 is a mobility, δE/δϕi are variational derivatives of the total energy E and are given by
The velocities of the components may be also determined in a thermodynamically and mechanically consistent way. Assuming that the solid and liquid volume fractions remain constant: and with constant in space and time, the resulting generalized Darcy laws for the velocities of the components are given by
where q is the water pressure, p is the solid pressure and 0, , j are positive definite motility matrices. The constitutive laws (59), (61) and (62) guarantee that in the absence of mass sources, the energy in equation (58) is non-increasing in time as the fields evolve. This approach may be straightforwardly extended to account for elastic and viscoelastic response of the solid fraction.
In related work Armstrong et al , and later Gerisch and Chaplain , considered a nonlocal model of adhesion in which j = 0, p = 0 and the variational derivatives in the velocity uj are replaced by (in 2D)
where n is the normal vector to a ball of radius R centred at x which is termed the sensing region and R is the sensing radius. The functions O and g determine the strengths of the interactions among the different cell types (and extracellular matrix) and ϕ = (ϕ1, . . . ϕN) is a vector of volume fractions. As the sensing radius R tends to zero, the nonlocal model converges to a local reaction–diffusion-taxis system of partial differential equations [39, 233]. For a finite sensing radius, Hillen et al  recently proved the global existence of solutions to the chemotaxis equations (see also the recent review by Hillen and Painter ). This adhesion velocity may actually be derived using energy variation arguments starting with a fully nonlocal version of the energy E such as
where Jij is an appropriately defined interaction kernel. See the appendix of  for a description of the procedure.
In other related work, Khain and Sander  developed a model of cell–cell adhesion for a single cell species similar to that described in Liquid–liquid model I above. In particular, Khain and Sander use a generalized Cahn–Hilliard  (GCH) equation to study tumour invasion. Below a critical level of adhesion, the tumour invades as a propagating front. Above the critical level, a second peak is found to appear in the tumour fraction (density) behind the leading front. The results of the GCH model compare well with a stochastic discrete cell model recently developed by Khain et al .
The model given in equations (59)–(62) may be simplified  by assuming (i) that tumour cells prefer to adhere to one another rather than to the host, as observed experimentally  and (ii) that no distinction is made between the adhesive properties of the viable and dead cells. Accordingly, in equation (53), we may take , where is the solid fraction of the tumour tissue and ϕN = ϕH is the volume fraction of the host tissue . Further, taking εij = 0 for i, j < N and , the total adhesion energy (58) is
This form of the energy arises also in the classic theory of phase transitions (e.g. ). For example, f may be written as the difference of the two convex functions
where one may take
where α1 and α2 describe the strength of adhesion (attraction) of tumour cells to the host tissue and each other, respectively, and Ē is an overall energy scale. Setting α1 = α2, yields a double well energy f with minima at and ϕT = 0 and gives rise to a well-delineated phase separation of the tumour () and the host tissues (ϕT ≈ 0). Since ϕT is continuous, it is necessary that in the interfacial region dividing the tumour and host domains. On the other hand, the states or ϕT < 0 are not physical, and the interaction energy tends to prevent their formation by increasing the energy of those states. Note that taking an interaction energy with logarithmic terms would explicitly prevent their formation . By varying the ratio α1/α2, the relative tendency of the tumour cells to aggregate can be modified. For example, increasing the ratio results in an increasingly diffuse tumour mass.
The thickness of the diffuse interface between the tumour and host tissue depends on the relative sizes of , Ē and α1/α2. Specifically, for fixed Ē the smaller the constant is, the less diffuse the interfacial region is. If Ē ~ 1/ε, then it can be shown that this system converges as ε → 0 to a classical sharp-interface single-phase tumour model . If the tumour contains different species that have different adhesion properties, the energy (65) can be modified to account for the different cell–cell interactions following the more general approach described earlier (see also ).
From the flux constitutive equation (59) and the adhesion energy (65), the adhesion fluxes may be determined. Recalling that the densities of the components are matched, and taking the mobilities i = ϕi, where is a positive constant, the fluxes are obtained :
where i = 1, . . . , N – 1 and , where it is used that the energy does not depend expicility on ϕH. The variational derivative is given by
Setting i = 0 for i > 0, which is consistent with assuming the cells are tightly packed and that they move together with the mass-averaged velocity, the component velocities become 
where it is used that the energy does not explicitly depend on ϕ0 and ϕN. In these equations, the terms dependent on δE/δϕT represent the excess force due to adhesion and arise from cell–cell interactions. The coefficients W, and are motilities that reflect the response of the water and cells, respectively, to pressure gradients. These coefficients may depend on the volume fractions and other variables as the individual components may respond to the pressure and adhesive forces differently, but mixtures of components tend to move together. The cell motilities contain the combined effects of cell–cell and cell–matrix adhesion. The constitutive choices (68), (71) and (70) guarantee that in the absence of mass sources (Si = 0), the adhesion energy is non-increasing as the fields evolve, while the total tumour mass is conserved.
As a first approximation, viable tumour cells may be assumed to necrose based only on the level of the local nutrient concentration n, i.e. when the nutrient n falls below the cell viability limit N,i which may be different for different cell types. Cells may be assumed to be comprised entirely of water. In terms of volume fraction, this is a reasonable first approximation. Cell mitosis may be assumed proportional to the amount of nutrient present and as mitosis occurs, an appropriate amount of water is converted into cell mass. Conversely, the lysing of cells represents a mass sink as cellular membranes are degraded and the mass converts completely into water. Mitosis is neglected in the host domain as the proliferation rate for tumour cells is much larger. Accordingly, defining ϕD to be the volume fraction of unlysed dead tumour cells, we may take 
where and are the rates of volume gain or loss due to cellular mitosis, apoptosis, necrosis and lysing, respectively, and where ∞ is the far-field nutrient level. Finally, is the Heaviside function. For simplicity, we have omitted mass exchange terms whereby mutations or epigenetic events may transform one cell species into another.
In the following discussion, ‘nutrient’ denotes a generic cell substrate, such as oxygen or glucose. Following the single-phase approach discussed in section 2.2, the host tissue is modelled at equilibrium as a first approximation, where the net nutrient uptake therein is negligible compared with the uptake by tumour cells. Nutrient uptaken by the host tissue is assumed to be replaced by supply from the normal vasculature. This may not be the case in the tumour, where not only the uptake exceeds the supply, but the uptake may also be much higher than that of the host tissue [179, 441]. Using the quasi-steady approximation described earlier for single-component tumour models, the nutrient transport equation may be written as
where D is the nutrient diffusivity and can vary depending on the medium and cell types, e.g. it can be an interpolated function from the host to the tumour, with constant value 1 inside the tumour and DH in the host medium. TC is a nutrient source from a pre-existing uniform vasculature or newly formed vasculature through angiogenesis and are uptake rates by the different cell types.
We now consider the growth of a tumour comprising of viable cells ϕV and dead cells ϕD, i.e. ϕT = ϕV + ϕD. The dead cell population includes the tumour cells that have undergone apoptosis or necrosis. Dead cells are assumed not to consume nutrient. There is also a host cell species and a liquid component. Using the generalized Darcy's law mixture model described in the previous section, the growth of an initially perturbed spherical tumour in 3D with cell adhesion γ = 0 is shown in figure 13. The isosurface ϕT = 0.5 is plotted. The parameters are the same as for figure 5 in . Substrate gradients contribute only to variable tumour cell proliferation and necrosis but do not induce migration (e.g. chemotaxis of tumour cells). Consistent with linear analyses [129, 339], the tumour is unstable due to the limited supply of nutrient and the corresponding tumour mass loss due to necrosis in the tumour interior. The simulation shows that in the nonlinear regime, a morphologically complex tumour emerges with repeated budding and sub-spheroid growth. By acquiring this morphology, the tumour has effectively increased its surface area to volume ratio in order to gain access to more nutrients from the surrounding host vasculature.
Cristini et al  considered the effects of chemotaxis induced through molecular species. Tumour evolution in nutrient-rich and nutrient-poor tissues was investigated by performing two-dimensional nonlinear numerical simulations. It was demonstrated that a tumour may suffer from taxis-driven fingering instabilities that are most dramatic when cell proliferation is low, as predicted by linear stability theory. This is also observed in experiments. This work shows that taxis may play a role in tumour invasion and that when nutrient plays the role of a chemoattractant, the diffusional instability is exacerbated by nutrient gradients. This model is thus capable of describing some of the complex invasive patterns observed in experiments .
Briefly, in Cristini et al  nutrient-driven taxis was included by adding the term χnnϕT to the total energy in equation (65), following the general fomulation given in equation (53). Equation (69) then becomes:
Using the liquid–liquid mixture model I, Cristini et al simulated the growth of a solid tumour in a nutrient-poor microenvironment. The proliferation rate is low and cell apoptosis is neglected. In figure 14(a), the tumour evolution from Cristini et al for low cell proliferation (we note that if χn = 0, i.e. when there is no taxis, the evolution is stable). The solid curve represents the tumour interface (ϕT = 0.5) and dash–dotted curve represents the linear results for the sharp-interface model. The figure illustrates that fingers develop around time t = 10 and get stretched out at time t = 20, forming long, slim and invasive fingers, thereby increasing the surface area of the tumour and allowing better access to nutrient. At later times, the fingers continue to stretch and the fingers tend to bend inwards. There is good agreement between the linear and nonlinear results at early times although there is significant deviation at later times due to strong nonlinearity. Figure 14(b) shows similar simulations of tumour evolution for a larger proliferation rate. All other parameters are the same as figure 14(a). The fingers are thicker for the case with larger proliferation and the spread of the fingers into the surrounding tissue for smaller proliferation is more pronounced at early times. As before, there is good agreement between the linear and nonlinear results at early times before nonlinear effects dominate the evolution.
The numerical results reveal that the tumour may suffer from nutrient-taxis-driven fingering instabilities which are most dramatic when cell proliferation is low. This suggests that nutrient-taxis may play a role in tumour invasion and that the diffusional instability is exacerbated by nutrient gradients. This resembles the branched invasive structures observed experimentally in tumour spheroids grown in hypoxic conditions .
The most common form of human glioma (cancer of the glial cells in the brain) is glioblastoma multiforme (GBM), which is the most aggressive stage of the disease. Typical survival rates of patients diagnosed with GBM are less than twelve months. GBM are characterized as being highly invasive, vascularized tumours (unlike some lower grade gliomas). See, for example, Preusser et al .
The results from the mixture model described in section 2.5.3, coupled with a lattice-free model of angiogenesis [417, 418] (see section 2.4.2), were recently compared with clinical brain tumour data by Frieboes et al , after estimating parameter values from in vitro  and in vivo experiments . Parameters governing the extent of neovascularization and nutrient supply due to blood flow were estimated in part from dynamic contrast enhanced magnetic resonance imaging (DCE-MRI) observations in patients  (see  for details). The model correctly predicted gross tumour morphology, including interior regions of necrosis surrounded by viable tumour tissue, as well as a tortuous neovasculature. Figures 15 and and1616 show a simulation of a growing glioblastoma using a recently updated version of this model (see ). Because the simulation was based upon the hypothesized relationships between cell adhesion, motility, death and proliferation described in , it allows for a quantitative analysis of these relationships in virtual patients, which could drive the development of morphological and immunohistopathological criteria for in vivo tumour invasion.
The model predicted regions of viable cells, necrosis in inner tumour areas and a tortuous neovasculature as observed in vivo . The vessels migrate towards the tumour/host interface since perinecrotic tumour cells and host tissue cells close to the tumour boundary produce angiogenic factors and other regulators. The tumour eventually coopts and engulfs the vessels.
The model also enabled prediction of the fine details of tumour morphology by quantifying cell metabolic response to spatial gradients in substrate levels caused by spatio-temporal heterogeneity in proliferation and the abnormal neovasculature. Excellent agreement was found when comparing the tumour ‘virtual histopathology’ with clinical pathology (figure 17), observing 100–200 μm thick sections of tumour encircling neo-vessels, surrounded by necrotic tumour cells and vessels that were shut down due to either age or mechanical intratumoural pressure (observed in ). This vessel shutdown enhances the hypoxic gradients seen clinically and predicted through computer modelling. The overall tumour shape depended strongly upon the vascular patterning, a result strongly supported by animal [57, 323, 327, 411, 459] and clinical imaging (e.g. ) and ex vivo histological data. The virtual histopathology also suggests that a tumour may rely upon vessels in the nearby host tissue [49, 426] and proliferate towards them; this is confirmed in the brain histopathology in figure 17.
The model enables a quantitative analysis, e.g., viable region thickness of about 100–200 μm and extent of necrosis as seen in figure 17 are shown to be strongly dependent on diffusion gradients of oxygen/nutrient in the microenvironment and agree with previous experiments [194, 269]. Further, the model predicts that the tumour boundary moves at a rate of about 50–100 μm per week, presenting a mass of diameter of about 5 cm in one year. These results are supported by well-known clinical observations (e.g. ). As the tumour grows and engulfs vessels in its vicinity, the tumour may compress the vessels  and disrupt flow of nutrients, leading to further necrosis and even temporary mass and vascular regression [275, 569]. In addition, chaotic angiogenesis leads to heterogeneous perfusion in the tumour that also may be responsible for regression of parts of the vascular network and necrosis of tumour cells [101, 407]. This enhances variable tumour cell proliferation.
The morphological instability leading to invasion in glioma resembles the overall shape of microsatellites observed experimentally in vitro and the predicted simulation morphologies in avascular conditions figure 18 (labelled ‘sub-spheroids’), including a smooth leading edge consisting of proliferating tumour cells and a trailing edge of necrotic cells. These images and data are in good agreement with morphologies predicted by stability analyses and computer simulations [126, 193, 194, 484].
These results support the idea that sophisticated multiphase tumour simulators, capable of simulating vascularized growth in two and three dimensions, and calibrated by in vitro and in vivo data, have the potential to predict tumour behaviour in patients.
Proliferation of epithelial cells that have undergone malignant transformation but remain at their original site, confined by their basement membrane, is called carcinoma ‘in situ’ . In particular, ductal carcinoma ‘in situ’ (DCIS) is the first stage of breast cancer, in which malignant cells proliferate inside the milk duct or lobule, increasing pressure on the basement membrane; at some point, they may breach it to invade the surrounding tissue . As tumour cells proliferate within the duct, they will eventually die. Calcium deposits form within the remains of these cells, which can be identified through imaging. Since mammography images may be difficult to interpret, Gavaghan et al  used breast carcinoma to illustrate how mathematical modelling can aid in detection, diagnosis and treatment. The authors discussed models underpinning mammography image analysis, which complement tumour growth models. They gave an overview of the primary image enhancement technologies, as well as a detailed description of their recent work in physics-based modelling in mammography. The goal of this theoretical approach to image analysis was to yield information that could be incorporated into the mathematical models.
Xu ( used the radially symmetric tumour model of Byrne and Chaplain  to study spatial patterns (e.g. stripes, spots and uniform distributions) shown by the stationary model solutions, noting that they were consistent with morphologies commonly observed in DCIS. In the model, diffusion of cell substrates limits tumour growth and the duct wall is assumed rigid; local pressure and cellular density are neglected. Franks and co-workers developed models that account for these effects [186–188], coupling existing models of avascular tumour growth in a cylindrically symmetric tube with mechanical models for the finite deformation of a compliant basement membrane. The coupling was mediated by interactions between the expansive forces of cell proliferation and the stresses that develop in the membrane (Epstein and Johnston ). Cell movement was described by a Stokes flow constitutive relation, and the effects of the material properties (i.e. the viscosity) on the tumour shape, as well as the extent to which cells adhere to the duct wall, were studied Franks et al . It was shown that stable, non-planar, interface configurations result; during the initial progression before the duct wall is breached, few cells die and a nutrient-rich model can be sufficient to capture the behaviour. In further work, Franks et al , interactions between the expansive forces created by cell proliferation and the stresses that develop in the compliant basement membrane were investigated, showing how the duct wall deforms during tumour growth, and how tumour progression along the duct depends on wall stiffness. Key model parameters were varied to determine how treatment, protease production and the inclusion of the surrounding stroma affect this growth.
These models provided insight into DCIS progression and generated hypotheses to be tested experimentally, e.g. the pressure on the duct wall is likely to be greatest at the tumour centre, indicating the likely location where the wall will be breached . The model was also used to test hypotheses for the localization of proteases that may compromise the wall. Model simulations and asymptotic analyses suggested that elevated pressure, rather than hypoxia, is more likely to stimulate protease production and localization near the wall . There is experimental evidence that in addition to stimulating protease production, local pressure may lead to basement membrane degradation through the death of myoepithelial cells, which are between the epithelial cells and the membrane .
The role of the microenvironment in breast cancer progression has been extensively studied by Gatenby and co-workers through both mathematical models and experiments (e.g. [181, 224–226, 228, 230, 239, 240). Although an evolutionary sequence of genetic malformation is observed in this progression, tumour growth is constrained by acidosis and hypoxia, which develop as cells proliferate into the avascular duct lumen. In order to survive, cells adapt to these adverse conditions by evolving acid-resistant and glycolytic phenotypes, which may be critical for invasion. The models show severe acidosis and hypoxia in regions more than 100 μm from the basement membrane of intraductal tumours, with invasive proliferation following the development of acid-resistant and glycolytic phenotypes. This hypothesis is supported by experimental and clinical evidence showing that growth into normoxic regions follows adaptation to hypoxia (e.g. ). An example of this is shown in figure 19 using a discrete cellular automaton model (described below).
In recent years, there has been much research using discrete methods to model tumour growth, invasion and angiogenesis. In this approach, individual cells or subcell elements are explicitly tracked and updated according to a set of biophysical rules. Discrete, or individual-based, models are generally divided into two categories: lattice-based (cellular automata) and lattice-free. In lattice-based modelling, the cells or subcellular elements are confined to a regular lattice. In lattice-free models, these elements may be placed anywhere throughout space. Typically, discrete models involve a composite discrete-continuum approach in the sense that substrate concentrations (e.g. oxygen, glucose, matrix-degrading enzymes, etc) are described using continuum fields while the cell-based components are discrete. See the reviews by Alber et al , Moreira and Deutsch , Drasdo , Araujo and McElwain , Quaranta et al , Hatzikirou et al , Nagy , Abbott et al  Byrne et al , Fasano et al , Galle et al , Drasdo and Höhme , Thorne et al , Anderson et al , Deisboeck et al , Anderson and Quaranta , Quaranta et al  and Zhang et al . Here, we briefly present several approaches to discrete modelling. Other examples and a more comprehensive list of references can be found in the review papers listed above.
In recent work, Rejniak [445, 446] developed a highly detailed approach to modelling the evolution of solid tumours. Each individual tumour cell is modelled using the immersed boundary method (e.g. see the review Peskin . The cell is represented as the interior of an elastic membrane with the nucleus represented as an interior point. In addition to the elastic forces, cell–cell adhesion and cell contractile forces are modelled using linear springs to mimic a discrete set of membrane receptors, adhesion molecules and the effect of the cytoskeleton on cell division, respectively (see figure 20). The cytoplasm and extracellular liquid are modelled as viscous, incompressible fluids. Elastic, adhesion and contractive forces impart singular stresses on the intra- and extra- cellular fluids. Cell proliferation is modelled by introducing a point source in the interior of the cell to increase its volume. Contractile forces act on opposite sides of the cell to create a neck that pinches off to produce two approximately equalsized daughter ‘drops’ (figure 21). Nutrient is supplied via diffusion and is modelled using continuum reaction–diffusion equations where the uptake term depends on the cell locations. In figure 22, the progression of an avascular tumour is shown. Note the proliferating cells (grey) at the (irregular) leading edge, followed by quiescent (white) and necrotic cells (black) in the interior. The model can describe individual cell morphology, although the method is computationally expensive thus restricting the number of cells to about 100. This work was later extended by Rejniak and Dillon  to incorporate a more realistic description of the cell membrane (two closed curves connected by springs) that better represents the lipid bilayer structure. In addition, sources and sinks are placed in the cell membrane to provide a model of water channels.
The cell-immersed boundary model has subsequently been used by Rejniak and Dillon to study preinvasive intraductal tumours, and by Rejniak and Anderson [447, 448], to study the formation and stability of epithelial acini. These structures form a basic building block of ducts in the breast and prostate, for example. Epithelial acini consist of a single layer of polarized endothelial cells enclosing a hollow lumen surrounded by a basement membrane and extracellular matrix. Genetic abnormalities in the endothelial cells may result in abnormal acini and ductal carcinoma. The immersed boundary model was further used by Dillon et al  and Anderson et al  to study the effects of nutrient availability, cell metabolism and phenotypic characteristics affecting the growth and stability of growing avascular tumours. The instabilities observed are in qualitative agreement with earlier continuum models (Cristini et al , Zheng et al , Li et al , Macklin and Lowengrub ) and cellular automaton models (e.g., Anderson , Anderson et al  and Gerlee and Anderson [234, 235]) discussed below.
A less detailed, but still cell-based approach, has been developed using the extended Q-Potts model which is adapted from Ising models in statistical physics and was originally developed to model surface-diffusion grain growth in materials science (e.g. Anderson et al ). Graner and Glazier [243, 247] adapted the Q-Potts model to simulate cell-sorting through differential cell–cell adhesion. In this approach, now referred to as the GGH (Graner–Glazier–Hogeweg) model, each cell is treated individually and occupies a finite set of grid points within a Cartesian lattice; space is divided into distinct cellular and extracellular regions. Each cell has a finite volume and a deformable shape. Cell–cell adhesion is incorporated through an energy functional. A Monte Carlo algorithm is used to update each Cartesian lattice point and hence change the effective shape and position of a cell. Although the description of the cell shape is less detailed than in the immersed boundary approach described above, finite-size cell effects are incorporated.
Stott et al  incorporated nutrient-dependent mitosis and necrosis in the GGH model to simulate the growth of benign, multicellular avascular tumours to a steady state. A continuum reaction–diffusion equation for nutrient was used. The steady-state configuration consists of an area of central necrosis surrounded by quiescent cells and an outermost shell of proliferative cells. Parameters are determined such that the thickness of these regions matches experimental results. Turner and Sherratt  later extended the GGH model by incorporating the ECM, the secretion of matrix-degrading enzymes, haptotaxis and adhesion-controlled mitosis to study tumour invasion. More generally, the GGH model has also been extended to account for chemotaxis , cell-differentiation  and cell-polarity . Jiang et al  further modified the GGH model to include a subcellular model of a protein expression regulatory network for the cell-cycle as well as incorporating growth promoters and inhibitors through continuum reaction–diffusion equations. Parameter values were determined such that the model produced avascular tumours that quantitatively replicate experimental measurements of spheroids in culture. The GGH model has also been used to simulate vasculogenesis by Merks et al [373–375], by including a model of contact inhibition, and by Bauer et al  in the context of tumour-induced angiogenesis in a heterogeneous tumour microenviroment. Rubenstein and Kaufman  recently studied the role of ECM in glioma invasion using a version of the GGH model. Scianna et al  used the GGH model to simulate cell scatter in MLP-29 mouse embryo liver cell tumours in response to hepatocyte growth factor. Very recently, Poplawski et al  studied the morphological evolution of 2D avascular tumours using the GGH model and developed a phase diagram characterizing the development of instabilities at the tumour/host interface and tumour morphologies with critical parameters measuring diffusional limitations (ratio of growth rate to diffusion rate) and cell–cell adhesiveness. In particular, they find that the development of instability depends primarily on the diffusional limitation parameter while the morphology details depend on cell–cell adhesion. The results are consistent with previous simulation results by Cristini et al [126, 129], Rejniak , Li et al , Macklin and Lowengrub , Gerlee and Anderson [234, 235] and Anderson et al .
Another Monte Carlo based model that incorporates finite cell size was developed by Drasdo et al  in the context of epithelial cell, fibroblast and fibrocyte aggregates in connective tissue. In this approach, the cells are simplified and consist of a roughly spherical space about a centre region. The cells are slightly compressible and are capable of migration, growth and division. An undividing cell is taken to be spherical. As a cell undergoes the mitosis process, it deforms into a dumbbell shape until its volume increases roughly by a factor of 2 and then divides into two daughter cells. Adhesion and repulsion (from limitations on cell deformation and compressibility) among cells are modelled using an interaction energy that describes nearest-neighbour interactions. Mitosis or cell migration may induce pressure on neighbouring cells. The cells respond by changing their mass or orientation to minimize the total interaction energy via a stochastic Metropolis algorithm . Note that in related work, Palsson and Othmer  modelled cells as deformable viscoelastic ellipsoids (see also Dallon and Othmer ).
Drasdo and Höhme  adapted the above Monte Carlo approach to avascular tumour spheroids in the early stages where growth is not primarily determined by oxygen or nutrient supply but is affected by volume exclusion due to limited cell compressibility. By comparing the results with experiments on tumour spheroids by Freyer and Sutherland  estimates were made for the biomechanical and kinetic parameters in the model. In later work, Drasdo and Höhme  extended the model to account for nutrient (glucose) limitations to growth, cell-necrosis and lysis of necrotic cells and simulated the spatio-temporal growth dynamics of two-dimensional tumour monolayers and three-dimensional tumour spheroids. Values of the model biophysical and kinetic parameters were drawn from the experimental literature. The model suggests that the transition from exponential to sub-exponential growth that is observed in experiments when the tumour is sufficiently large is due to biomechanical growth inhibition. Glucose deprivation was found to primarily determine the size of the necrotic core but not the size of the tumour. Galle et al  extended the model to incorporate cell–substrate contact-dependent cell death (anoikis ) and studied monolayer cell growth of epithelial cells. They found that inactivation of cell–substrate contact-dependent cell-cycle arrest, cell–substrate-dependent programmed cell death or cell–cell contact mediated growth inhibition can lead to epithelial tumour growth. In the context of growing monolayers, Drasdo  showed how the finite-sized Monte Carlo approach can be used to determine the rules for a simpler cellular automaton model (in which cells have zero size, see below) and how the cellular automaton model can be used to derive a continuum model with contact inhibition by coarse-graining thereby providing a link between different scales and biophysical processes. A similar study was performed by Byrne and Drasdo  who presented further analysis of the continuum model. Radszuweit et al  used these models to analyse the effect of geometry and investigated the differences in growth in cell populations in 2D and 3D. Höhme and Drasdo  investigated the effects of biomechanics and nutrient to determine their relative importance on the growth of cell populations. Ramis-Conde et al  incorporated the dynamics of E-cadherin–β-catenin interactions to obtain a more realistic model of cell–cell adhesion (E-cadherins are cell–cell adhesion proteins; β-catenin binds the membrane-bound E-cadherins to the cell cytoskeleton). This detailed model is capable of describing detachment of cells from a primary tumour and the corresponding epithelial-to-mesenchymal transition. Ramis-Conde et al  extended this model to simulate tumour cell intravasation into blood vessels and investigated the role of cadherins in metastasis. Galle et al  incorporated cell–matrix interactions and stroma contact-dependent cell regulation as well as cell-differentiation. They studied the effect of cancer stem cell organization on tumour growth. They found that when stem cells are located on the tumour/host boundary this may result in more rapid invasion of the host tissue than if the stem cells were confined to the tumour interior.
Following earlier work by Honda et al , Meineke et al  and Brodland and Veldhuis , Schaller and Meyer-Hermann [470, 471] modified the approach described above and represented cells as Voronoi polyhedra in three dimensions (e.g. polygonal boundaries) instead of deformed spheres thus enabling a more realistic description of densely packed cells. Schaller and Meyer-Hermann also incorporated a description of the cell cycle into their model and considered the effects of both oxygen and glucose. By matching their results with experiments on tumour spheroids, they find that the ratio of oxygen to glucose uptake ratio should be about 1 : 5; their model predicts the number of cells in each phase of the cell cycle for different concentrations of oxygen and glucose.
Models in which cells are represented by discrete points (i.e. variations in cell size and shape are not considered) give a still less detailed discrete approach. Cells can be interpreted as having a fixed size (e.g. ranging from 5–20 μm in radius). These methods can be both lattice-based (cellular automata) and lattice-free (agent-based). An advantage of these algorithms is that due to their simplicity, they are easier to implement and can be solved more efficiently than the more detailed discrete methods described above. Here, we present a few recent examples. Further references may be found in the review papers listed above.
Kansal et al [303, 305] developed a cellular automaton model to simulate the growth of 2D and 3D brain tumours using an adaptive Delaunay triangulation (Voronoi tessellation) to minimize grid-induced anisotropy. The density of lattice sites varies throughout space (e.g. larger site density in the centre of a tumour) where the lattice size reflects the number of cells represented. Simulations were performed for avascular, multicell spheroids and reproduce a Gompertzian growth pattern. The authors then used this approach to investigate the emergence of subclonal subpopulations . In Gevertz and Torquato , this approach was extended to incorporate a dynamic model of angiogenesis where a random analogue of a Krough cylinder model was used to model the regression and sprouting of a tumour-induced neovasculature. The vasculature development is coupled with the growth of the tumour. Simulations were performed of vascular brain tumour growth at early stages. This work was later extended by Gevertz et al  to simulate tumour growth in confined, heterogeneous environments.
Anderson et al  developed a cellular automaton model of solid tumour growth focused on four variables related to tumour invasion: tumour cells, host tissue (ECM), matrix-degrading enzymes and oxygen. The latter three are continuous (concentrations) while tumour cells are discrete (individuals). Cells move via a biased random walk on a Cartesian lattice; the coefficients of a discretized version of the partial differential equation governing the tumour cell density to generate probabilities of movement. The transition probabilities are similar in spirit to the gradient model developed earlier by Othmer and Stevens  in the context of cell-chemotaxis. The cells produce matrix-degrading enzymes and respond haptotactically to the density of extracellular matrix. Cell–cell adhesion was not considered. The model predicts a greater extent of local tumour invasion in a heterogeneous extracellular matrix than is predicted by the analogous continuum model. This approach has been generalized by Anderson and coworkers to include cell–cell adhesion, different clones (mutations) and phenotypes as well as other effects and processes (see below).
Patel et al  described a cellular automaton model to simulate early tumour growth, and examined the roles of host tissue vascular density and tumour metabolism in the ability of a small number of monoclonal transformed cells to develop into an invasive tumour . The model was two-dimensional and incorporated continuum hydrogen ion and glucose fields obeying partial differential equations. This work was later extended by Smallbone et al ) and Gatenby et al  to incorporate different phenotypes and mutations. Smallbone et al show that space and nutrient limitations are sufficient to drive the emergence of hyperplastic, glycolytic and acid-resistant phenotypes that have the capability to breach the basement membrane and form an invasive tumour.
Schmitz et al  studied the effects of treatment by extending an automaton model of brain tumour growth, finding that spatial location, i.e. confinement, downgraded the competition of resistant clones with respect to the sensitive cells. As a result, the fraction of resistant cells is a less important indicator of patient prognosis when the cells are confined. If mutations arise in response to treatment, tumours with both very frequent and very infrequent mutations develop with more spherical geometries, while tumours with an intermediate level of mutations display multi-lobed geometries, where mutant strains are at localized points on the surface of the tumours. Mansury et al  incorporated cell migration and a simple model of mechanical resistance to movement in a cellular automaton model of tumour growth. By allowing cells to search both locally and globally for a more growth-promoting environment, the model shows that a phase transition may occur such that if the search is global, a few large clusters form with a short lifespan while if the search is local, many small slower-growing and longer-lived clusters form.
Dormann and Deutsch  developed a lattice-gas cellular automaton method for simulating the growth and size-saturation of avascular multicell tumour spheroids (see also the book  by the same authors). Unlike traditional cellular automaton methods where one cell can be at a single grid point (volume exclusion), lattice-gas models accommodate variable cell densities by allowing multiple cells at a single grid point by allowing separate channels of movement. The channels specify direction and velocity magnitude, which may include zero velocity resting states. Lattice-gas models typically require channel-exclusion-only one cell may occupy a channel at any time. This approach has later been used by Hatzikirou et al  to study travelling fronts characterizing glioma cell invasion. Very recently, Hatzikirou et al  performed a mean field analysis of a lattice-gas model of moving and proliferating cells to demonstrate that certain macroscopic information (e.g. scaling laws) can be accurately obtained from the microscopic model. Accurate predictions of other macroscopic quantities that sensitively depend on higher order correlations are more difficult to obtain.
Chen et al  developed a stochastic model for the growth of a heterogeneous solid tumour by deriving a master equation that accounts for two tumour-cell types. Explicit expressions governing the variances and covariance were obtained and simulations were performed to determine the evolution of the cell populations and their uncertainty.
Alarcón et al  developed a mathematical model to show how blood flow and red blood cell heterogeneity influence growth of normal and cancerous cells; these were considered as elements of a cellular automaton, with evolution rules determined from different behaviour of normal and cancer cells. These authors also used a physiologically structured lattice model for vascularized tumour growth that accounts for interaction between cancerous and normal tissue, cell division, apoptosis, blood flow and structural adaptation of the vasculature, transport of oxygen and release of vascular endothelial growth factor . Effects of nutrient heterogeneity as well as growth and invasion of cancerous tissue were investigated. This work was later extended by Owen et al  to model tumour growth coupled with dynamic network models for the vasculature.
Building on the earlier work of Anderson et al (2000), Anderson (2005) extended the cellular automaton model to include cell–cell adhesion. Cell–cell adhesion is modelled by weighing the probability of motion by the number of desired neighbours. Different cell adhesion phenotypes are incorporated by different numbers of desired neighbours and proliferation and nutrient uptake rates. The microenvironment plays an important role through cell–extracellular matrix (ECM) interactions. In the model, cells degrade the ECM which produces and maintains nutrient gradients through uptake by cells. This model enabled evaluation of how individual cell– and cell–matrix interactions may affect the tumour shape. This work was extended  to provide a theoretical/experimental framework to quantitatively characterize tumour invasion as a function of microenvironmental selective factors. Mutations are assumed to occur randomly; in , 100 different phenotypes are considered. Cell genotypes may evolve to support invasion driven by definable microenvironmental selective forces. In agreement with the findings of Cristini et al [126, 129, 194, 575], conditions such as hypoxia and heterogeneous extracellular matrix were found to induce invasiveness through fingering tumour margins, and dominated by more aggressive cell phenotypes. The evolution of an avascular tumour using this approach is given in figure 23. The colour gradiations correspond to different cell phenotypes (necrotic–black, dark grey—highly proliferative and invasive (low adhesion), lighter grey—quiescent). Gerlee and Anderson  used a simpler version of this model to investigate complex branched growth patterns that arise during cell colony growth under nutrient limited conditions. In agreement with earlier stability analyses (e.g. ) the stability of the growth was found to depend on how far the nutrient penetrates into the colony. For low nutrient consumption rates the penetration distance was large, which stabilized the growth, while for high consumption rates the penetration distance was small, which led to unstable branched growth. By extending this approach to include a feed-forward neural network to model the decision-making mechanisms governing the evolution of cell phenotype, Gerlee and Anderson [234, 236] demonstrate how the oxygen concentration may significantly affect the selection pressure, cell population diversity and morphology of a tumour. Gerlee and Anderson  later extended this model to study the emergence of a glycolytic phenotype. Their results suggest that this phenotype is most likely to arise in poorly oxygenated tissue with large matrix density. More recently, Piotrowska and Angus  extended this model to incorporate multiple cells at each cellular automaton lattice site and studied the effects of necrosis in avascular tumour growth. In particular, their model suggests that acidity is not sufficient to initiate necrosis before nutrient levels drop below levels needed for cell viability.
Bartha and Rieger  developed a cellular automaton algorithm to track coupled tumour growth and tumour-induced neovascularization using a discrete approach for both. Each lattice point may be occupied by a tumour cell and/or a blood vessel with a particular radius which is determined by the blood flow rate and wall shear stress. The evolution of each follows a stochastic process involving proliferation, migration and death on a rectangular grid. Nutrient and growth factor concentrations are modelled using approximations of Green's functions to the diffusion operators. Simulations reveal that vascular collapse in a few critical segments can dramatically influence the flow through the network and the tumour progression. This work was used by Lee et al  to develop a theoretical model of tumour vascularization using percolation theory. They suggested that vascular collapse in the tumour centre leads to a percolation process driven towards criticality by vessel stabilization due to increased blood flow in the remaining vessels. Welter et al  extended the model by incorporating the Fahraeus effect and refining the simulation of oxygen. Welter et al  found that vessel regression and collapse may increase perfusion due to the fact that their tumour was vascularized primarily co-option of the existing vasculature. Welter et al  recently extended the model to account for an arterio-venous network.
Mantzaris  developed a Monte Carlo algorithm to simulate the dynamics of heterogeneous cell populations, taking into account the intrinsic randomness of cell division and the unequal partitioning of cellular material at cell division. The effects of genetic mutations in a tumour cell population have also been studied, for example, as a continuous model of cell densities coupled with a discrete model of the cell cycle . Mansury and Deisboeck  presented a numerical agent-based model of tumour cells to study the molecular, microscopic and multicellular patterns that emerge from various interactions among cells and their environments. Mallet and Pettet  developed a model of haptotaxis which includes a description of integrin activation, local recruitment and protrusion from the cell body. In the model, the cell aggregate was then able to respond to a gradient of cell–matrix adhesion, represented by functionally active integrins.
Wang et al  developed a multiscale model to study the growth of non-small cell lung cancer within a 2D microenvironment, implementing a specific intracellular signal transduction pathway between the epidermal growth factor receptor (EGFR) and extracellular receptor kinase (ERK) at the molecular level. Phenotypic changes at the cellular level were triggered through dynamical alterations of these molecules. The results indicated that for this type of cancer, downstream EGFR–ERK signalling may be processed more efficiently in the presence of a strong extrinsic chemotactic stimulus, leading to a migration-dominant cell phenotype and an accelerated rate of tumour expansion. Zhang et al  presented a three-dimensional multiscale agent-based model to simulate the cellular decision process to either proliferate or migrate in the context of brain tumours. Each cell was equipped with an EGFR gene–protein interaction network module that also connected to a simplified cell-cycle description. The results show that proliferative and migratory cell populations directly impact the spatio-temporal expansion patterns of the cancer. This was later refined by Zhang et al  to incorporate mutations representing a simplified tumour progression pathway. Very recently, Quaranta et al  examined interactions between the tumour microenvironment and cancer cells at various scales using cell-based discrete mathematical models.
Ramis-Conde et al  developed an individual-based lattice-free model of solid tumour growth where the tumour cells interact with one another through a potential function which incorporates cell-attraction and repulsion. The cells respond by moving down gradients of the potential but may also undergo haptotaxis up gradients of extracellular matrix and chemotaxis up gradients of a chemical agent, both of which are modelled as continuum variables. Simulations show finger-like invasions of tumour cells where the leading cells degrade the extracellular matrix to form a channel and the trailing cells follow a gradient of chemoattractant released by the degraded ECM in their wake. Another stochastic agent-based model of tumour growth was developed recently by Gomez-Mourelo et al . Their model extends the Anderson cellular automaton model by removing the lattice-based constraint, and incorporating an age-structure. Gomez-Mourelo et al derive a macroscopic system that describes the limiting behaviour of the discrete model, when the number of cells tends to infinity. Simulations show that that the limiting macroscopic system of partial differential equations and the agent-based model yield similar results for tumour fingering patterns.
An important research direction for the future involves the development of hybrid continuum–discrete models for tumour growth. Hybrid models have the potential to combine the best features of both continuum and discrete approaches, and may provide more realistic coupling of biophysical processes across a wide range of length and time scales. As mentioned earlier, there are two types of continuum–discrete models. In the first, which we term composite models, tumours are represented as the sum total of discrete elements (e.g. cells) while cell substrates (e.g. nutrients, growth and other chemical factors) and ECM are represented as continuum variables. There have been many models developed using this approach (see section 3). In the second type, which we term hybrid models, the tumour itself is described using both continuum and discrete components. There have been far fewer efforts to simulate tumour growth using the second approach. However, we believe that hybrid modelling is very important as it affords the possibility of seamlessly upscaling from the cell-scale to the tumour and tissue scales. Because of the computational cost, it would be very difficult to reach this range of scales using models where the tumour cell colonies are only described using discrete cells although there are now techniques for accelerating these methods such as the heterogeneous multiscale method (e.g. ) and the equation-free method (see below).
In an early study, Othmer and Stevens  derived general classes of partial differential equations to study the continuum limit of reinforced random walks of discrete cells where the cells move in response to chemical signals, and found that a variety of dynamics could occur in a cell population, e.g. aggregation, blowup or collapse. As a way to bridge multiple time scales, Setayeshgar et al  developed a numerical evolution scheme for a class of stochastic problems in which the temporal evolution occurs on widely separated time scales, and for which the slow evolution can be described in terms of a small number of moments of an underlying probability distribution. Setyehsgar et al (2005) applied this approach to simulate the motion of non-interacting bacteria cells in the presence of a chemoattractant. This work combined a kinetic Monte Carlo method, to model bacteria motility, and a continuum ordinary differential equation inner solver, to model signal transduction internal to the bacteria, with a projective integration outer solver. The authors showed that projective time integration of coarse variables (e.g. bacteria densities) could be carried out on time scales long compared with that of the microscopic dynamics. In particular, numerical results and analysis of errors in support of the efficacy of this integration were presented. The algorithm provided a means to integrate the effective macroscopic equations, which were not known explicitly.
Connecting the cell-scale with the tumour scale is a challenge in cancer modelling. Erban and Othmer  provided a connection between the micro- and macro-scales describing bacterial chemotaxis by deriving the macroscopic description from a microscopic model of the behaviour of individual cells. In the model, each bacterium was assigned an internal state that evolves according to a system of ordinary differential equations forced by a time-and/or space-dependent external signal. The turning rate was a function of the internal state of the cell, which in turn depends on the external signal. Solutions of the macroscopic equations agreed well with Monte Carlo simulations of individual cell movement. Erban et al  then used simplified versions of these continuum and discrete models to test the efficiency and accuracy of equation-free methods to accelerate the discrete random walks models. The equation-free approach, developed earlier by Kevrekidis et al  (see also the review by Li et al ) leverages spatio-temporal scale separation to allow significant gains in computational efficiency by alternating short, localized bursts of appropriately initialized microscopic simulations with accelerated processing of the results (e.g. temporal ‘jumping’) at the macroscopic, continuum scale. In particular, certain macroscopic variables such as the cell density may vary smoothly over relatively long space and time scales while macroscopic details have rapid fluctations over small space and time scales. Although the macroscopic equations for the cell densities are known in the specific example of chemotaxing cells considered by Erban et al (2006), the equation-free approach may be applied to systems in which the coarse-scaled continuum governing equations are unavailable and too difficult to explicitly derive. Erban et al (2003) obtained a decrease in computation time by over a factor of 1000 using the equation-free approach to simulate large number of chemotaxing cells. Later, Bold et al  proposed an equation-free approach to model the collective dynamics of populations of coupled, heterogeneous biological oscillators.
Recently, Kim et al  and Stolarska et al  coupled an agent-based algorithm, to model cells in the outer proliferating rim of an avascular tumour, with a continuum description of cells in the inner quiescent and necrotic regions of the tumour. The extracellular matrix and cell substrates (e.g. oxygen and glucose) are also modelled as continuum fields. The idea is that a millimetre-sized tumour spheroid contains on the order of 106 cells, which is too many to feasibly simulate using a direct individual-based approach. On the other hand, the proliferating rim of cells is much smaller—on the order of 100–200 μm in thickness—and contains about 10 times fewer cells. While this is still a significant number, a discrete simulation is somewhat more feasible if one uses an accelerating algorithm such as the equation-free approach described above. In the much larger quiescent and necrotic regions, it is much more efficient to use a continuum description.
Kim et al defined four geometrically distinct regions. These are, from the outside to the inside of a tumour, the extracellular matrix or culture medium surrounding the tumour, a region of proliferating cells at the outer edge, a quiescent zone bordering the proliferating region, and a necrotic core (see figure 24). The cell-based component described the mechanical interaction of cells with the surroundings, how an individual cell (modelled as an ellipsoid) reacts to forces on it, and how growth and division are described and the effect of stress on them. The cell-based component extends an earlier model developed by Dallon and Othmer  to include growth and cell division. The cell description accounts for forces the cells induce on one another and the extracellular matrix, dynamic drag forces that arise as a cell changes the number of attachments with neighbouring cells and a static frictional force that models rigid attachment to other cells or the extracellular matrix. Forces are transferred from the discrete cell system to the continuum by interpolation from the discrete cells to the nearby nodes of a mesh triangulation on which the continuum problems are posed. A linear viscoelastic model is assumed to hold for the stress in the continuum regions. As discrete cells in the proliferative region become quiescent, and vice-versa as quiescent cells become proliferative, a least-squares projection algorithm is used.
In figure 24 evolution of a tumour spheroid in the absence of an outer gel is shown from Kim et al (2007). The discrete cells, numbering 253 at t = 43 h, are indicated by the dark small circles; the spatial unit is 10 μm. Adjacent to the discrete cells is a quiescent region and necrosis occurs in the inner (white) region. As the tumour grows, slight asymmetries form because of variations introduced in part by the sequential updating of cell states in addition to the variations in nutrient fields due to slightly different cell sizes that drive asymmetries in proliferation, quiescence and necrosis rates. By treating the proliferating cells using a discrete model, this approach enables the use of greater biological detail than would be the case had the description of cell growth been fully continuous. For example, the effect of variations in cell-cycle time, intra and intercellular mechanics and adhesiveness, cell size and shape can be examined.
In very recent work, Bearer et al  and Frieboes et al (in preparation)  have developed an alternate approach to tumour modelling using both continuum and discrete descriptions of tumour cells. In their approach, both the continuum and discrete representations of tumour cells are used simultaneously throughout space, subject to mass and momentum conservation laws which incorporate interactions among the discrete and continuous fields. An agent-based, lattice-free model is used for discrete cells. The discrete cells are assumed to have zero size. A circular random walk model, that incorporates chemotaxis, haptotaxis and volume exclusion such that two cells cannot occupy the same position in space, is used for the discrete cells. The cells may also respond to the velocity induced by the continuum volume fractions and compete with the continuum description for nutrients. Accordingly, additional spatially localized uptake terms are introduced in the nutrient transport equation (76) to represent nutrient uptake by the discrete cells. The continuum volume tumour cell volume fractions are obtained from a mixture model of the type described in section 2.5.3 (see also Wise et al .
In Bearer et al and Frieboes et al (in preparation), the discrete cells arose in hypoxic regions of the tumour due to the downregulation of cell–cell and cell–matrix adhesiveness. This may represent an epithelial-to-mesenchymal transition from collective motion to individual motion and is consistent with experimental observations of the effects of hypoxia [197, 198, 267, 302]. Bearer et al did not account for any interaction forces in the discrete and continuum cell descriptions. However, as mentioned above, discrete cells may respond to the velocity induced by cell proliferation and cell adhesion in the continuum tumour fields.
The mass conservation equation for the continuum cells accounts for the exchange of mass from the continuum-to-discrete representations (and vice-versa). Assuming that there is a single clone of viable tumour cells as was done in Bearer et al (2009), mass is exchanged from the continuum to the discrete fields via a rate of exchange function. Suppose it is determined that at a location xi a discrete cell should be created. For example, in Bearer et al (2009), such locations are identified in regions of hypoxia. The rate of exchange of mass from the continuum to the discrete field may be taken to be 
where Cextract is a rate of extraction, ϕV is the volume fraction of viable tumour cells and G(x − xi) is a function localized around the region x ≈ xi. The total volume transferred from the continuous to the discrete field in a time t, measured from the time tdetect when it is determined that a discrete cell should be created, is . Let the time t = tcreate denote the time that Vc,i(tcreate) = Vcell, the target volume for a discrete cell. At this time, a discrete cell is then created at the location xi. An analogous mass-conserving algorithm is used to convert discrete cells into the continuum volume fraction field. In regions where the density of discrete cells exceeds a threshold value due to either increased proliferation or aggregation, the discrete cells are converted back to the continuum volume fraction. Accordingly, the following source function is added to the continuum tumour cell equation:
where Cdeposit is a rate of conversion and the sum is taken over all cells in the high density region. A cell is removed at time tremove when reaches the cell target volume: Vr,i(tremove) = Vcell. Both the continuum and discrete cells uptake nutrients and oxygen and thus an additional uptake term is incorporated to account for uptake by discrete cells.
In figure 25, the evolution of a vascularized tumour using the scheme described above is shown. Discrete cells (small blue dots) are released from hypoxic (perinecrotic) regions of the continuum tumour description (grey surfaces). The cells move away from the tumour bulk and up gradients of oxygen, which is supplied at the far-field boundary and from the newly formed vasculature (dark brown curves). Vessel sprouts are also shown (yellow curves). A lattice-free angiogenesis model is used and coupled to the tumour growth model following the methods described earlier by Frieboes et al  in the context of glioma modelling (section 2.5.5). As the discrete cells evolve, they degrade and remodel the extracellular matrix by laying down fibronectin macromolecule networks. The cells respond to the remodelled microenvironment by forming single-file like strands of palisading cells from low oxygen (perinecrotic regions) to high oxygen environments. As the cells reach oxygen-releasing vessels, their proliferation is upgraded and their motility is downgraded. As their density exceeds the threshold, the discrete cells are converted back to the continuum and tumour microsatellites form around the neovasculature. The vasculature also responds to the cell-proliferation pressure and vessels are shut off and the vasculature regresses if the cell-pressure is sufficiently large. This can be seen to occur around some of the newly formed continuum clusters. The resulting morphology compares well with pathology data  for human brain tumour specimens.
In the simulation described above, the discrete cells interact with each other directly only through a volume exclusion velocity and indirectly through the nutrient field as both continuum and discrete cells uptake nutrients (and the availability of nutrients limits proliferation). To incorporate further interactions between the continuum and discrete descriptions, we may derive interactions among the continuum and discrete forces by introducing an energy that accounts for both continuum and discrete effects. Then, by following the multiphase modelling approach outlined earlier, we may naturally derive thermodynamically consistent interaction forces. Consider the total (nonlocal) energy of the system to be
where Jij is an interaction potential and ϕi, ϕj are volume fractions of the i and j components, which may be either continuum or discrete (cf equation (64)). Let us suppose that ϕ0, . . . , ϕN represent continuum fields while ϕN+1, . . . , ϕM represent discrete fields. In particular, for the discrete fields ϕj(x) = wjδ(x − xj) where j > N, wj is a weight, δ is the delta function and xj denotes the location of the cell. Then, we may rewrite the energy as
where Ecc represents the self-interactions among the continuum fields
the energy Edd represents the self-interactions among the discrete cells
while the remaining energy describes the interactions between the discrete and continuous components
One may now derive constitutive relations for cell velocities and fluxes using an energy variation argument. Following an approach analogous to that used for the second liquid–liquid mixture model in section 2.5.2, we obtain exactly the same generalized Fick's law as in equation (59) for the continuum fluxes and the same generalized Darcy law velocities as in equation (62) for the continuum volume fractions, but with the total energy E given above. In particular, the variational derivatives contain the interaction terms:
where the former is fully continuous but the latter contains the effects of the discrete cells. For example, assuming for simplicity that Jij(x, y) = Jij(x − y) is an even function and Jij = Jji, we obtain
For the discrete cells, we may take the over-damped evolution
where ζi > 0 is a mobility and the overdot denotes the time derivative. This is similar in spirit to the lattice-free discrete model developed recently by Ramis-Conde et al  described earlier where cells interact through a discrete potential, see section 3. Here, the first term δEdd/δxi characterizes the discrete potential while the second contains the effects of the continuum:
With these choices, the total energy E is non-increasing in time in the absence of mass sources (thermodynamic consistency). The effects of haptotaxis and chemotaxis may be added (as in the simulation shown in figure 25) and these models may be simplified by idealizing the interactions among the different species. For example, as mentioned above, Jij may be highly localized and some interactions among the cell species may be neglected. The models may also be generalized to account for elastic and viscoelastic forces in the continuum description by considering the corresponding system energy, and by generalizing the constitutive relations and the continuum dynamics.
In this paper we have provided a limited overview of recent results on theoretical cancer modelling. As we have shown, there have been significant efforts to model cancer as a complex system where there is variability and coupling among biophysical processes across a wide range of spatial and temporal scales. While we have discussed the use of discrete models, we have focused more on the continuum approach, and we have reviewed recent examples of the incorporation of biologically relevant parameter values into multiscale models of tumour growth and invasion. We have also reviewed hybrid frameworks where the tumour tissue is modelled containing both discrete (cell-scale) and continuum (tumour-scale) elements, thus connecting the micrometre to the centimetre tumour scale. While computational models and simulation results still lag behind experimentation in studying cancer progression in vivo, significant advances have been made in the nonlinear modelling of the progression of this complex disease.
Mathematical modelling predicts that transport limitations of cell nutrients, oxygen and growth factors may result in cell death that leads to morphological instability. This can provide a mechanism for invasion via tumour fingering and fragmentation. These conditions induce selection pressure for cell survivability, and may lead to additional genetic and epigenetic responses. The results suggest that the tumour morphology and dynamics are coupled in complex, nonlinear ways to cell phenotype and to molecular properties (e.g. genetics) as well as to phenomena in the environment such as hypoxia. These properties and phenomena act both as regulators of morphology and as determinants of invasion potential by controlling cell proliferation and migration mechanisms [198, 483, 529]. The importance of this close connection between tumour morphology and the underlying cellular/molecular scale is that through mathematical modelling (e.g. morphology) may be used to understand the underlying cellular physiology and predict invasive behaviour and response to treatment.
Although we did not discuss models of cancer therapy here, these results have important implications for therapy since decreasing nutrient levels in the microenvironment, or inducing large-scale cell death, tends to increase tumour fragmentation and invasion into the surrounding tissue as well as drug resistance (e.g. [191, 465, 485]). Indeed, several experimental studies have recently shown that hypoxia, for example, as induced by anti-angiogenic therapies may result in the development of multifocal tumours [57, 151, 327, 454, 459, 462, 474]. By combining therapy treatments with anti-invasive drugs such as metastatic (Met) inhibitors [48, 69, 382] or hepatocyte growth factor (HGF) antagonists [145, 377] the development of multifocal tumours can be greatly reduced . This is confirmed by increasing the cell–cell adhesion (or decreasing cell mobility) in mathematical models (see figure 12 and ).
Conversely, increasing nutrient levels in the tumour microenvironment may lead to greater morphological stability of tumours making it easier to remove surgically. This suggests  that treatments that seek to normalize tumour vasculature (by selectively ‘pruning’ weak blood vessels with targeted anti-angiogenic therapy ) may stabilize tumour morphology by providing increased access to nutrient. Since such treatments may also increase the accessibility to chemotherapeutic agents [292, 484], mathematical analyses provide additional support for the use of targeted anti-angiogenic therapy as adjuvant to chemotherapy and resection.
While great strides have been made in recent years on the development of multiscale models of solid tumour growth, there is much more work to be done to provide a more comprehensive understanding of cellular diversity and adaptation by describing the complex interactions among tumour cell clones and their microenvironment [483, 529]. An important direction for solid tumours concerns the role of stem cells and their lineages, which were not discussed in this paper. Increasingly, it is being recognized that solid tumours are organized in a hierarchical manner through cell lineages arising from collections of cancer stem cells. In this picture, stem cells divide to self-renew and to produce a lineage of more differentiated cells with (possibly) limited replication potential. Thus stem cells are hypothesized to play critical roles in tumourigenesis, tumour progression and recurrence after treatment. See the recent reviews by Pardal et al , Wicha et al , Dalerba et al , and Visvader and Lindeman . While there have been numerous mathematical models of cancer stem cells in liquid (hematopoietic) tumours such as leukemia, there are far fewer models of stem cells in solid tumours as evidence for the existence of stem cells in solid tumours has only been recently established. In the past several years, mathematical models incorporating stem cells have been developed for colorectal cancer, e.g. see van Leeuwen et al [530, 531] and Johnston et al [296, 295], as well as for breast cancer, e.g., see Enderling et al  and Ashkenazi et al . Generic tumour models incorporating cancer stem cells have also been developed, e.g. see Ganguly and Puri , Piotrowska  and Galle et al . Endering et al  used a cellular automaton model to investigate the effects of migration, proliferation and death rates on cell lineages in cancer, finding that migration of stem cells is a key variable that enhances the growth of a tumour via the spatio-temporal seeding of new colonies (termed self-metastasis). This model was also used to study the effect of the response of tumour colonies to radiotherapy .
For further references on mathematical modelling of cancer stem cells we refer the reader to the recent review by Michor . However, many questions remain regarding the effect of spatio-temporal heterogeneity in the microenvironment and the response of the cells to regulating chemical signals (e.g. see Jackson et al  for a discussion of the role of cell signalling models in cancer). The links among molecular, cellular and tissue-scale processes need to be further explored and integrated with experimental observations and clinical data.
Realistically, linking across such a wide range of scales requires the use of hybrid continuum–discrete methods and holds considerable promise for the future. These methods will allow the possibility of further expanding the current understanding of invasion and migration [124, 170, 195, 196, 198, 307, 319, 452, 461, 483, 529, 556, 565] and instead describe the progression as the result of a complex system of nonlinearly interacting multiscale processes. It is with this in mind that the prediction of disease progression and treatment response through mathematical modelling based on patient-specific tumour characteristics can become a reality.
The authors thank the anonymous referees for their insightful comments. JSL, SMW and VC gratefully acknowledge partial support from the National Science Foundation (NSF) Division of Mathematical Sciences. JSL and XL are also grateful for partial support from the National Institutes of Health (NIH) grants NIH-1RC2CA148493-01 and NIH-P50GM76516 for a National Center of Excellence in Systems Biology at UCI. VC also acknowledges partial funding from Cullen Trust for Health Care, Department of Defense, and the NIH.