PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Nonlinearity. Author manuscript; available in PMC 2010 August 30.
Published in final edited form as:
Nonlinearity. 2010; 23(1): R1–R9.
PMCID: PMC2929802
NIHMSID: NIHMS223837

Nonlinear modelling of cancer: bridging the gap between cells and tumours

Abstract

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.

1. Introduction

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 [246], 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 [198] 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 [492] 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 [326], 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 [255]; (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 [6], Chaplain [106], Bellomo and Preziosi (2000), Moreira and Deutsch [381], Bellomo et al [59], Swanson et al [507], Araujo and McElwain [34], Mantzaris et al [362], Friedman [199], Ribba et al [450], Quaranta et al [439], Hatzikirou et al [260], Nagy [389], Byrne et al [88], Fasano et al [182], Jackson et al [286], van Leeuwen et al [531], Roose et al [456], Graziano and Preziosi [248], Harpold et al [256], Friedman et al [201], Sanga et al [464], Martins et al [365], Deisboeck et al [150], Anderson and Quaranta [30], Bellomo et al [60], Cristini et al [127], Wang and Deisboeck [541], Preziosi and Tosin [430], Jackson et al [534] and Tracqui [525]). 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 [317], Stolarska et al [498] and Bearer et al [55]).

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.

2. Continuum modelling

2.1. Background

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, 508511, 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, 8487, 92, 93, 95, 96, 120, 129, 146, 186, 189, 220, 222, 272, 283285, 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 ([117119, 198]). ECM reorganization by tumour cells [198] 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. [457]), viscoelastic (e.g. [346, 317], and elasto-viscoplastic (e.g. [24]. Theoretical nonlinear analyses of various tumour models mentioned above have been performed (e.g. [53, 78, 122, 131134, 136142, 155, 182, 199, 202204, 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.

2.2. Basic tumour model

2.2.1. Overview

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 [482] accounted for a mitotic inhibitor and analysed the stability of growth, McElwain and Morris [368] accounted for apoptosis, and Adam [6] 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 [105] presented mathematical models of spherical tumour growth through the stages of avascular growth, angiogenesis and vascularization, as well as pattern formation in cancer [106]. Friedman and Reitich [206, 207] and Cui and Friedman [138] 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 [129] 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 [129] thus suggesting that the mechanisms that control tumour morphology also control its ability to invade. Recently, Li et al [339] 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]).

2.2.2. Mathematical formulation

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 [251], Byrne and Chaplain [89], Friedman and Reitich [206], Cristini et al [129], 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 [nabla] · u is given by

u=λp,
(1)

where λp is the cell-proliferation rate and λp is given by

λp=bnλA,
(2)

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.

Figure 1
Diagram of a basic non-necrotic tumour. The tumour occupies the volume Ω, Σ is the interface between tumour tissue and health tissue, n is the unit outward normal to Σ.

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:

0=D2n+Γ,
(3)

where Γ is the rate at which nutrient is added to Ω and is given by

Γ=λB(nBn)λn,
(4)

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.

To determine the cell velocity, Darcy's law may be used as the constitutive assumption (e.g. [89, 129, 206, 250]):

u=μP,
(5)

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. [575]). Models of viscoelasticity (e.g. [346]), elasto-viscoplasticity [16] 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

(n)Σ=n,
(6)
(P)Σ=γκ,
(7)

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

V=μn(P)Σ.
(8)

Following [89, 129, 206, 250] and others, assume that λ, λA, λB, nB, b are uniform. Following [129], denote λM = bn to be the characteristic mitosis rate, λR=μγLD3 to be the intrinsic relaxation time scale, and B = nBλB/nB + λ) to be a measure of the extent of vascularization. Introducing the non-dimensional length scale LD=D12(λB+λ)12, and time scale λR1, a modified concentration Γ and pressure [p with macron] can be defined [129]:

n=n(1(1B)(1Γ)),P=γLD(p+(1Γ)G+AGxx2d),
(9)

where G and A measure the relative strength of cell–cell and cell–matrix adhesion and apoptosis, respectively:

G=λMλR(1B),A=λAλMB1B.
(10)

Dropping the bar notation, the non-dimensional equations for Γ and p can be obtained:

2ΓΓ=0,
(11)
2p=0,
(12)

with boundary conditions:

(Γ)Σ=1,
(13)
(p)Σ=κAG(xx)Σ2d
(14)

in a d-dimensional tumour (d = 2, 3). The non-dimensional normal velocity of the tumour–host interface is

V=n(p)Σ+Gn(Γ)ΣAGn(x)Σd.
(15)

2.2.3. Regimes of growth

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

dRdt=V=AGRd+G{I1(R)I0(R),d=2,(1tanh(R)1R),d=3.}
(16)

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.

(1) Low vascularization

G ≥ 0 and A > 0 (i.e. B < λAM). 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 [251] 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 [91].

(2) Moderate vascularization

G ≥ 0 and A ≤ 0 (i.e. 1 > B ≥ λAM). 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 VG as R → ∞.

(3) High vascularization

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: λAM > 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).

2.2.4. Linear analysis

Consider a perturbation of the radially symmetric tumour interface Σ:

rΣ=R(t)+δ(t){cos(lθ),d=2,Yl,m(θ,ϕ),d=3,}
(17)

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]:

(δR)1d(δR)dt={Gl(l21)R3Gl+2RI1(R)I0(R)GIl+1(R)Il(R)I1(R)I0(R)+lAG2,d=2,Gl(l+2)(l1)R3(Gl+3R+GIl+32(R)ll+12(R))(1tanh(R)1R)+lAG3,d=3.}
(18)

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).

Linear stationary states

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 [129] 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.

Linear self-similar evolution

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,

A(l,G,R)={2(l21)GR3+2(l+2l)1RI1(R)I0(R)+2lIl+1(R)I1(R)Il(R)I0(R)2l,d=2,3(l1)(l+2)GR3+3(l+3l)1R(1tanhR1R)+3lIl+32(R)Il+12(R)(1tanhR1R)3l,d=3.}
(19)

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.

Figure 2
Analysis of the basic tumour. Top: apoptosis parameter A as a function of unperturbed radius R from condition (19) for self-similar evolution; d = 2 (dashed) and d = 3 (solid); G and l labelled. Asymptotic behaviours are dotted (see [129]). The two solid ...

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.

Diffusional instability

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 [385]. In the high-vascularization regime, shrinkage may be unstable.

2.2.5. Model calibration

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 [194], 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 [384]. 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 = (Doxyoxy,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 [102] and Doxy = 1.45 × 10−5 cm2 s−1 [393] 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 PLD2λM3μ 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 [194].

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 [194] 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 [129]. 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.

Figure 3
Evaluation of tumour stability. Diagram shows death parameter A versus spheroid radius R (rescaled with diffusion length L). The curves for given values of G are obtained from [194] governing spheroid morphological stability. Experimental conditions for ...

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.

2.2.6. Nonlinear results

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 [206] and Cui and Friedman [139] for radially symmetric solutions using the Darcy law model; Cui [131] for a model with inhibitors; Friedman and Reitich [207], Bazaliy and Friedman [53, 54], Cui [132], Cui and Wei [141] and Cui [134] for non-symmetric solutions; Friedman [200] and Wu and Cui [558] for a Stokes model; Zhou, Escher and Cui [577] for a model of multilayer tumours and Cui [135] for a hyperbolic equation model of tumour growth. Bifurcations from spherically symmetric solutions have been studied by Friedman and Reitich [208], Friedman and Hu [203205] and Zhou and Cui [576]. In recent work, Xu and Gilbert [564] 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 [129] and Li et al [339]) 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.

Nonlinear stationary states

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 [129], 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 G=GlNL<Gl, where Gl is the result from linear theory (section 2.2.4). It is found that GlNLGl=O(δ2), where δ is a measure of the perturbation size. Thus, nonlinearity is destabilizing for the stationary shapes.

Nonlinear self-similar evolution

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: R=Areaπ. 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.

Figure 4
Analysis of the basic tumour. Left: self-similar shrinkage for R0 = 4 and δ0 = 0.2 (t = 0 to 0.96 shown). Right: Unstable shrinkage for R0 = 4 and δ0 = 0.4 (t = 0 to 0.99). The solid curves correspond to nonlinear solution and the dashed ...

Nonlinear evolution in the high-vascularization regime

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. [105]). 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 [394] of in vivo angiogenesis and tumour growth.

Nonlinear unstable growth in the low vascularization regime

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 [129]. 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).

Figure 5
Evolution of the basic tumour surface in the low-vascularization regime, for d = 2, A = 0.5, G = 20 and initial tumour surface. Dotted lines: solution from linear analysis; solid: solution from a nonlinear calculation with time step Δt = 10−3 ...

An analogous evolution is observed in 3D [339]. 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.

Figure 6
Evolution of the basic tumour surface in the low-vascularization regime, A = 0.5, G = 20 and initial tumour surface as defined in [129]. (a) t = 0, (b) t = 2.21, (c) t = 2.42, (d) t = 2.668. Reprinted with permission from Li et al Discrete Contin. Dyn. ...

The biological significance of the diffusional instability [129] 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]).

2.3. Tumour growth in heterogeneous tissue

2.3.1. Overview

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 [561]. In addition, the microenvironment may induce clonal diversity [308].

2.3.2. Hypoxic microenvironments

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 [91] 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 [222] 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 [partial differential]Ω, 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

0=(Dn)λn(x,n)n+λbulkn(1n)B(x)χH,
(20)

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 [352]. These domains may actually be defined in terms of the oxygen concentration as in equation (21). In equation (20), λbulkn 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

λn(x,n)={0ifxΩH1ifxΩP={x1n>np}λQnifxΩQ={xnpn>nN}λNnifxΩN={xnNn},}
(21)

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 λNn. 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

(μP)={0ifxΩHnAifxΩP0ifxΩQGNifxΩN,}
(22)

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 [129] and others, cell–cell adhesion forces are modelled in the tumour by generalizing the condition (7) and introducing a Laplace-Young jump condition:

[P]=G1κ,xΣ,
(23)

where κ is the mean curvature and G=λMλR. 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.

PP,x(ΩTΩH).
(24)

2.3.3. Nonlinear results

In Macklin and Lowengrub [351], an accurate ghost-cell [183, 244]/level-set [397] 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 [198], 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 [350]. 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.

Figure 7
Tumour morphological response to the microenvironment. The external tissue nutrient diffusivity D increases from left to right and the external tissue mobility μ increases from bottom to top. The shape of each tumour is plotted at time T = 20.0. ...

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 [129] 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 [126] who also found that anti-invasive therapy (increased adhesion, decreased mobility) may be used successfully as adjuvant to anti-angiogenic therapy.

Complex tissue

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 [351]). 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.

Figure 8
Tumour growth in heterogeneous tissue. Simulation from time t = 0.0 days (top left) to t = 60.0 days (bottom right) in 10 day increments. White band on the right of each frame models a rigid material such as the skull; black denotes an incompressible ...

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 [347350]) 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.

Extracellular matrix, taxis and invasion

The effects of extracellular matrix on tumour invasion were modelled by Anderson et al [29] 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 [254] 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 [369] showed that in chemotaxis of cells in symmetric multicell spheroids may act against cell-motion induced by pressure gradients. Castro et al [103] 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 [314] 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 [311] find that symmetric fronts become unstable leading to instabilities and fingering as observed in experiments (e.g. [149, 152]). Chaplain and Lolas [112] developed a continuum model of ECM degradation by a particular matrix degrading enzyme (urokinase plasminogen) and studied the resulting tumour progression. Marchant et al [363] 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 [130], Woodward et al [557], Burgess et al [80], Swanson et al [508512] and Jbabdi et al [293]). 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 [512] 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 [387] that describe interactions between the tumour and host cells, Gatenby et al [229] 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 [578] proved the existence of travelling wave solutions to a system of equations originally proposed by Sherratt [479] and used later by Chaplain and Sherratt [480]. In this approach, nonlinear diffusion is introduced to model contact inhibition between tumour and host tissues.

Boushaba et al [71] 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 [39] 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 [520] analysed a mathematical model of avascular tumour spheroid growth that accounted for cell-cycle dynamics and chemotaxis-driven cell movement. Gerisch and Chaplain [233] 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 [513]. In [352], 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

u=μ(E)P+χEE.
(25)

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 [14] analysed cellular traction using an inverse problem approach. This was later used by Ambrosi et al [16] 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 [117119] 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 [514] 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. [297]). 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 [180]. 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).

Acidosis

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 [228] and Gillies et al [241]. 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 [224] 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:

Lt=rNVdL+DL2L,
(26)

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 [224] 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 [226].

Smallbone et al [488] 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 [226] 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).

Clonal diversity

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 [389]. Gatenby and Frieden [223] demonstrated this using information theory and showed that this is required in order to prevent the accumulation of genetic malfunctions. Gatenby and Vincent [231] 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 [240] 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.

Internal stress

Experiments demonstrate that mechanical stress may significantly influence the growth of solid tumours (e.g. [268]). The mathematical modelling of residual stress development in growing tumours was recently reviewed by Araujo and McElwain [34], Roose et al [456] and Tracqui [525]. Briefly, the study of Shannon and Rubinsky [477] 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 [114] to classify solid tumours and by Chen et al [120] 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. [268]).

Jones et al [298] 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

0=2nn,
(27)
v=nA,
(28)
σ=0,
(29)
12(v+vT)=13(v)I+12(DDt(3σTr(σ)I)+3(ωσσω)),
(30)

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 ω=12(vvT) and I the identity tensor. The Jaumann derivative has been used to maintain frame-indifference (D/Dt = [partial differential]t + v · [nabla] 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 [457]), by incorporating anisotropy (see Araujo and McElwain [35, 38]) or viscoelasticity (e.g. MacArthur and Please [346]). The effect of solid stress on the vasculature was investigated by Sarntinoranont et al [468] using the poroelastic model and by Araujo and McElwain [36] 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 [24] point out, the strain is not frame invariant [356] and thus the Jaumann derivative is inappropriate. Instead, the use of accretive forces (e.g., [13, 1922, 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 [24], the above models (with the convection terms dropped) may be obtained from a small deformation limit of the constitutive laws derived in [24] 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:

F=FrG,
(31)

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 [67] 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 [22] derived the following equations for a tumour spheroid

ρ=ρ0Jr1,
(32)
(F1(Jn(F1)T))=(nA)ρJ,
(33)
P=0,
(34)
g.=g3(nA),
(35)

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 = u, 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 [22] 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 [19], 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 [156] 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 [13] 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 [343] 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.

2.4. Tumour growth and neovascularization

2.4.1. Background

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 [255]. 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 [101]. 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.

Figure 9
Schematic (not to scale) of a necrotic tumour in transition from avascular to vascular growth. Disjoint regions ΩH, ΩV and ΩN represent healthy tissue viable tumoural tissue and necrotic core domains, respectively. Tumour region ...

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 [258], which leads to increased flow resistance and ultimately to heterogeneous supply of cell substrates in the tumour microenvironment [209].

2.4.2. Modelling of angiogenesis

Mathematical models of tumour-induced angiogenesis date to the work of Balding and McElwain [47]. 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 [90], Orme and Chaplain [396], De Angelis and Preziosi [146], Sansone et al [467], Levine et al [334, 335, 337], Hogea et al [272], Peterson et al [413], Stamper et al [491], Jain et al [288] and Addison-Smith et al [9]. 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 [497], Anderson and Chaplain [28], Tong and Yuan [522], Plank and Sleeman [417, 418], Sun et al [502, 503], Kevrekidis et al [310], Bauer et al [51], Milde et al [379] and Capasso and Morale [100]. Blood flow and network remodelling have also been simulated (e.g. Pries et al [437], McDougall et al [367], Stephanou et al [494, 495], McDougall et al [366], Wu et al [559], Zhao et al [574], Sun and Munn [501] and Pries and Secomb [435]). In addition, models of tumour growth in static network topologies have been simulated (e.g. Alarcón et al [10], Betteridge et al [66]) and multiscale models of fluid transport in tumour have also been developed (Chapman et al [116]).

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 [575]. This was modelled in 2D using a version of the continuum–discrete model of Anderson and Chaplain [28] coupled to the non-symmetric tumour growth model of Cristini et al [129]. 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 [238], Bartha and Rieger [50], Bauer et al [51] and Scislo and Dzwinel [547]. The effects of blood flow through the neovascular network on tumour growth were also recently considered in Alarcón et al [11], Bartha and Rieger (2006), Lee et al [330], Welter [548] and Owen et al [399] 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 [549]. Very recently, Macklin et al [352], 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 [366] to explicitly analyse the effects of blood flow and vessel remodelling. In other recent work, Lloyd et al [342] 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 [343]). In 3D, vascular tumour growth has been studied using a mixture model and lattice-free description of tumour-induced angiogenesis (Frieboes et al [193] and Bearer et al [55]).

Basic angiogenesis model

A basic angiogenesis model can be constructed [575] based on that of Anderson and Chaplain [28] with some additions taken from Paweletz and Knieriem [409], Paku [402], Chaplain and Stuart ( [115]) and McDougall et al [367]. The model of Anderson and Chaplain [28] 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’ [242], 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: [28]

  • concentration of tumour angiogenic factor (TAF, e.g. vascular endothelial cell growth factor VEGF) c
  • concentration of matrix-degrading enzymes (MDE) m
  • endothelial cell density (ECD) e
  • density of ECM (e.g. matrix macromolecule fibronectin) f.

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 [352]

0=(DTT)λdecayTT(λbindingTT)1sprout tips+λprod.T,
(36)

where DT=DT is the diffusion coefficient, 1sprout tips is the characteristic function of the sprout tips, λprod.T,λdecayT,λbindingT denote the non-dimensional production, natural decay and binding rates of TAF, and λprodT 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 [partial differential]T/[partial differential]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 [28]

ft=ηpeηmfm,
(37)

where ηp is rate of production of fibronectin by ECs and ηm is the rate of degradation of fibronectin by MDEs [575]. The MDE concentration, m, satisfies

mt=DmΔm+λprod.m(1m)1ΩV+λspr.prod.m1sprout tipsλdecaymm,
(38)

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 [575] 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 [82].

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 [575]:

et=De2e((χC1+αcc+χff+χuu)e)ρDe+ρpe(1e)H(cc)ρNχΩNe
(39)

where De is the EC diffusion coefficient, χc, χf are chemotaxis and haptotaxis coefficients, respectively, χu and α are dimensionless coefficients and H 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 [28], 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 [366] 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 [575]

dxdt=μCu,
(40)

where x is the position on the capillary and μC is the capillary mobility [575]. 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 [419] 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 [310] and Whitaker et al [551] incorporated the effects of inhibitors, following earlier work by Levine [337], in continuum–discrete models of angiogenesis. Bauer et al [51] modelled angiogenesis using the Graner–Glazier–Hogeweg (GGH) framework of cell-level modelling (see section 3). Milde et al [379] incorporated the effects of both matrix-bound and soluble growth factors in a three-dimensional model of angiogenesis. More recently, Capasso and Morale [100] 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 [436], Pries et al [432, 433]. Following this work, McDougall et al [367] 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 [494] 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 [495] (due to shear and circumferential stresses generated by flowing blood [433, 434]). This model was further updated by McDougall et al [366] 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 [495] and Alarcón et al [10]).

Coupling tumour growth with angiogenesis

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 [575] and Macklin et al [352], 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 λneoσ is added to the right-hand side of equation (20)

λneoσ=λneoσBneo(x,t)H(hHDhmin)(1C(Pvessel,P))(1n),
(41)

where λneoσ 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 [H with macron]D and hmin 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 [366]. This provides one aspect of the coupling between the tumour growth and angiogenesis models. A second mode of coupling between the two models occurs through the cutoff function C (Pvessel, P) 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 λneoσ, 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 [11].

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 [575] and Macklin et al [352], 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 [352] the source term λprod.T in equation (36) describes the source of TAF due to secretion by hypoxic cells in the perinecrotic region:

λprod.T=λprod.T(1T)1ΩQ,
(42)

where 1ΩQ is the characteristic function of the region of quiescent (hypoxic) cells and λprod.T denotes the non-dimensional production rate of TAF.

2.4.3. Simulation of neovascularization, diffusional instability and tissue invasion

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 [352]. 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.

Figure 10
Tumour-induced angiogenesis and vascular tumour growth. The vessels respond to the solid pressure generated by the growing tumour. Accordingly, strong oxygen gradients are present that result in strongly heterogeneous tumour cell proliferation and shape ...
Figure 11
Dimensional intravascular radius (m) and pressure (Pa) along with the non-dimensional ECM and TAF concentrations from the simulation shown in figure 10. Reprinted with permission from Macklin and Lowengrub J. Math. Biol. 58 788. Copyright © 2009 ...

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 [348351].

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 [126] 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’ [126] (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 [198] 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 [492]. 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.

Figure 12
Effects of treatment on tumour morphological stability. Solid red: calculated tumour boundary, black: necrosis. The neovasculature (pink) forms from ‘free’ endothelial cells (blue) that chemotax after sprouting from pre-existing vessels ...

Following the strategy of Zheng et al [575], Frieboes et al [193] simulated angiogenesis and vascular tumour growth in 3D using a multiphase tumour model [555] coupled with a lattice-free continuous–discrete model of angiogenesis originally developed by Plank and Sleeman [418]. 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.

2.5. Multiphase modelling

2.5.1. Overview

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 [246], 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 [73] 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 [533] and Bertuzzi et al [65] to account for ATP production from energy metabolism in the models of multicellular tumour spheroids.

Ambrosi and Preziosi [23] reviewed several mixture models and developed a mixture model treating the tumour as a deformable porous material. Byrne et al [95] and Byrne and Preziosi [97] extended the two-phase model of Ambrosi and Preziosi (2002). Jackson and Byrne [287] and Lubkin and Jackson [345] 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 [74] 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 [457] 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 [24]. Chaplain et al [110] 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 [43] developed a fully multiphase mixture model of tumour cords growing around a capillary. Later, Preziosi and Tosin [431] 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 [44] used a two population model to investigate the upregulation of glycolysis (Warburg effect) in tumour cords. Ambrosi and Preziosi [24] used a mixture model to study the elasto-viscoplastic response of the tumour and microenvironment during growth. Preziosi et al [428] 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 [249] 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 [216] 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 [34], Hatzikirou et al [260], Quaranta et al [439], Byrne et al [88], Graziano and Preziosi [248], Astanin and Preziosi [42], Roose et al [456], Preziosi and Tosin [430] and Tracqui [525]. 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 [555] developed a multispecies mixture model and performed simulations of tumour growth in three dimensions. Frieboes et al [193] and Bearer et al [55] 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 [98] 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 [39], and used later in the tumour context by Gerisch and Chaplain [233]. We also note that a continuum model of cell adhesion and migration was recently developed by Kuusela and Alt [325] 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, 553555]. The model was employed to study 3D vascularized cancer growth (Frieboes et al [192]), including malignant brain tumours (Frieboes et al [193], see also Sanga et al [464]).

2.5.2. General conservation equations

The primary dependent variables in a (N + 1)-species model can be specified as follows [555]:

  • the volume fractions of the water, tumour and host cell species, ϕ0, . . . , ϕN,
  • the densities of the components ρ0, . . . , ρN,
  • the stresses σ0, . . . , σN
  • the component velocities u0, . . . , uN.

It may be assumed that there are no voids (i.e. the mixture is saturated) and thus i=0Nϕi=1. 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

ρ(ϕit+(uiϕi))=Ji+Si,
(43)

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 u=i=0Nϕiui. Summing equation (43), the mass of the mixture is conserved only if i=0NJi is constant. Without loss of generality, this constant is chosen to be zero:

i=0NJi=i=0NSi=0.
(44)

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

0=σi+πi,
(45)

where πi are interaction forces among the species. Assuming that the mixture stress σ=i=0Nσi satisfies [nabla] · σ = 0, this gives the constraint i=0Nπi=0.

Let ui be the internal energy of the ith component and let the volume-averaged internal energy of the mixture be u=i=0Nϕiui. Then, in its simplest form, the balance of energy equation is given by (e.g. [37, 128])

ρϕiDiuiDt=σi:vi+ρiϕiri+βi,
(46)

where Di/Dt = ([partial differential]/[partial differential]t) + ui · [nabla] 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

i=0N(βi+Siuiπiui)=0.
(47)

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 η=i=0Nϕiηi. 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 [526]

ρ(ηt+uη)+Jρrθ0,
(48)

where J is an entropy flux and r=i=0Nϕiri, 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 [123].

A solid/liquid biphasic model

For binary mixtures of one solid (i = 1) cellular component and water (i = 0), Araujo and McElwain [37], assume that ψi = ψi(F1, G1, u1u0). In this approach, F1 is the deformation gradient in the solid phase and G1 = [nabla]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

σ0=ϕ0pI,σ1=ϕ1pI+ϕ1F1(ψ1F1)T,
(49)

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:

ϕ1ρu1=S1,ρD1ϕ1Dt=0(ϕ1u1+ϕ0u0)=0,
(50)

for conservation of mass. The growth equations are

DE1Dt=u1Ω+12μDDt(σ1+13trσ1)+ϕ12μDpDt(3ΩI),
(51)

where DDt is the Jaumann derivative (e.g. [300]), Ω is an anisotropic growth tensor [35] 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

S1=ϕ1ρDgDt3ϕ12ρ2μDpDt.

Finally, the momentum equations are

σ1+α(u0u1)=0,ϕ1p=α(u0u1),
(52)

where α is a drag coefficient.

Liquid–liquid mixture model I

An alternative approach was taken by Cristini et al [128]. 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

Ψi=Ψ~i(ϕ0,,ϕN)+ϕil=1Lχilcl+j=0Nil22ϕj2,
(53)

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 [555]). 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

ϕit+(ϕiui)=(SiJi)ρ,i=0N(ϕiui)=0,
(54)

for the conservation of mass. In equation (54), the adhesion fluxes may be given as [128] Ji = −M[nabla]μi which is a generalized Fick's law, for i ≠ 0, where μi is the chemical potential

μi=Fi(ϕ0,,ϕN)+l=1L(χilχ0l)clj=0N(jiΔϕij02Δϕ0),
(55)

and Fi=j=0N((Ψ~jϕi)(Ψ~0ϕ0)). The flux in the liquid is J0=i=1NJi. In the liquid (assumed to be inviscid), the momentum balance equation is the multicomponent Darcy's law:

ϕ0p=i=0Nαi(uiu0),
(56)

where the αi are drag coefficients. For the solid components, the general viscous law is obtained (for i > 0):

αi(uiu0)=ϕi(p+μi)+(Li(ui+uiT)),
(57)

where Li(ui+uiT)=λi(ui+uiT)+νi(ui)I, 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 = −i[nabla]μi, see [128]. Viscosity can be easily included and yields a nonlocal equation analogous to the mixture model derived in Byrne and Preziosi [97] in one-dimension; the inviscid model is similar to that considered earlier by DeAngelis and Preziosi [146]. 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 [128]. A discussion of the model without the above approximation may also be found in [128].

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[nabla]μi and that the velocity of the liquid is u0=1ϕ0i=1Nϕiui. This is analogous to the closure laws described in DeAngelis and Preziosi (2000) and Ambrosi and Preziosi [23].

Liquid–liquid mixture model II

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

E=i=0NΨidx,
(58)

where the energy density is still given by equation (53). As derived in [555], thermodynamically consistent fluxes may be taken to be the generalized Fick's law:

Ji=Mi(δEδϕiδEδϕN),1i<N1,
(59)

and JN=i=1N1Ji, where [M with macron]i > 0 is a mobility, δE/δϕi are variational derivatives of the total energy E and are given by

δEδϕi=j=0N(Ψ~jϕi+l=1Lχjlcl(εji2ϕi)),1i<N.
(60)

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: i=1Nϕi=ϕ~s and ϕ0=1ϕ~s with ϕ~s constant in space and time, the resulting generalized Darcy laws for the velocities of the components are given by

u0=k0(δEδϕ0+q),
(61)
uj=k(pi=1NδEδϕiϕi)kj(δEδϕj1ϕ~Si=1NϕiδEδϕi+pϕ~S),j1,
(62)

where q is the water pressure, p is the solid pressure and k0, k, kj 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 [39], and later Gerisch and Chaplain [233], considered a nonlocal model of adhesion in which kj = 0, p = 0 and the variational derivatives in the velocity uj are replaced by (in 2D)

A[ϕ](x,t)=1R0Rr02πn(θ)O(r)g(ϕ(t,x+rn(θ)))dθdr,
(63)

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 [271] recently proved the global existence of solutions to the chemotaxis equations (see also the recent review by Hillen and Painter [270]). This adhesion velocity may actually be derived using energy variation arguments starting with a fully nonlocal version of the energy E such as

E=12i,jJij(x,y)g(ϕi(x))g(ϕj(y))dxdy,
(64)

where Jij is an appropriately defined interaction kernel. See the appendix of [555] for a description of the procedure.

In other related work, Khain and Sander [312] 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 [98] (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 [313].

2.5.3. A special case of the liquid–liquid mixture model II

The model given in equations (59)(62) may be simplified [555] by assuming (i) that tumour cells prefer to adhere to one another rather than to the host, as observed experimentally [40] and (ii) that no distinction is made between the adhesive properties of the viable and dead cells. Accordingly, in equation (53), we may take Ψ~i(ϕ0,,ϕN)=ϕif(ϕT), where ϕT=i=1N1ϕi is the solid fraction of the tumour tissue and ϕN = ϕH is the volume fraction of the host tissue (ϕT+ϕH=ϕ~S). Further, taking εij = 0 for i, j < N and NN=ε, the total adhesion energy (58) is

E=Ω(f(ϕT)+ε22ϕT2)dx.
(65)

This form of the energy arises also in the classic theory of phase transitions (e.g. [98]). For example, f may be written as the difference of the two convex functions

f(ϕT)=fc(ϕTϕ~S)fe(ϕTϕ~S),
(66)

where one may take

fc(ϕ)=E4α1((ϕ12)4+1)andfe(ϕ)=E4α2(ϕ12)2,
(67)

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 ϕT=ϕ~S and ϕT = 0 and gives rise to a well-delineated phase separation of the tumour (ϕTϕ~S) and the host tissues (ϕT ≈ 0). Since ϕT is continuous, it is necessary that 0<ϕTϕ~S<1 in the interfacial region dividing the tumour and host domains. On the other hand, the states ϕT>ϕ~S 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 [169]. 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 [555]. 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 [316]).

Fluxes and velocities

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 [M with macron]i = [M with macron]ϕi, where [M with macron] is a positive constant, the fluxes are obtained [555]:

Ji=MϕiδEδϕT,
(68)

where i = 1, . . . , N – 1 and JN=JH=i=1N1Ji, where it is used that the energy does not depend expicility on ϕH. The variational derivative is given by

δEδϕT=f(ϕT)ε22ϕT.
(69)

Setting ki = 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 [555]

u0=uW=kWq,
(70)
ui=k(pδEδϕTϕT),
(71)

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 kW, and k 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.

Mass exchange terms

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 with macron]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 [555]

Si=λM,innϕiλA,iϕiλN,iH(nN,in)ϕi,N>i>0andiD,
(72)
SD=i>0iDN1(λA,i+λN,iH(nN,in))ϕiλLϕD,
(73)
SN=SH=0,
(74)
SW=i=1NSi=i=1iDN1λM,innϕi+λLϕD,
(75)

where λM,i,λA,i,λN,i and λL are the rates of volume gain or loss due to cellular mitosis, apoptosis, necrosis and lysing, respectively, and where [n with macron] is the far-field nutrient level. Finally, H is the Heaviside function. For simplicity, we have omitted mass exchange terms whereby mutations or epigenetic events may transform one cell species into another.

Nutrient diffusion

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

0=(D(ϕ1,,ϕN)n)+TCi=1iDN1νiuϕi,
(76)

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 νiu are uptake rates by the different cell types.

2.5.4. Nonlinear results

Morphological instability: simulations of a four-phase model

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 [555]. 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.

Figure 13
Multiphase model: growth of a two-plus-four mode tumour in 3D with cell adhesion γ = 0. The ϕV = 0.5 isosurface is shown. Parameters are the same as for figure 5 in [555]. The model predicts that this tumour morphological instability, ...

Chemotaxis-driven instability

Cristini et al [128] 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 [411].

Briefly, in Cristini et al [128] nutrient-driven taxis was included by adding the term χnT to the total energy in equation (65), following the general fomulation given in equation (53). Equation (69) then becomes:

δEδϕT=f(ϕT)ε22ϕT+εχnn.
(77)

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.

Figure 14
Effects of chemotaxis on tumour morphology: evolution of the tumour surface, with A=0,D=1, ε = 0.005, χσ = 5 and the initial tumour surface as (x(α) 12.8, y(α) − 12.8) = (2 + 0.1 cos 2α)(cos α ...

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 [411].

2.5.5. Comparison with clinical data

Glioma

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 [426].

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 [193], after estimating parameter values from in vitro [194] and in vivo experiments [193]. 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 [395] (see [193] 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 [192]). Because the simulation was based upon the hypothesized relationships between cell adhesion, motility, death and proliferation described in [193], 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.

Figure 15
Study of human glioma using a multiscale 3D mixture model [192, 193, 555]. Viable tissue (region between red outer tumour boundary and inner purple boundary), necrosis (region interior to inner purple boundary) and vasculature (thick red lines: mature ...
Figure 16
Study of human glioma using a multiscale 3D mixture model [192, 193, 555]. Another view of the simulation from figure 15 at a slightly later time (90 days).

The model predicted regions of viable cells, necrosis in inner tumour areas and a tortuous neovasculature as observed in vivo [79]. 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 [400]). 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. [560]) 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.

Figure 17
Study of human glioma using a multiscale 3D mixture model [192, 193, 555]. (a): Details of virtual tumour histology showing invasive tumour front (white); locations of blood-conducting new vessels (NV) and non-conducting vessel sprouts (blue dots). Aged ...

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. [388]). As the tumour grows and engulfs vessels in its vicinity, the tumour may compress the vessels [400] 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].

Figure 18
Simulations and experiments of cell protrusions in glioma growing into detached cell clusters and forming separate tumours. Top row: simulation snapshots (time =days) Bottom row: in vitro. Bar: 130 μm. Tumour microsatellites are seen on right. ...

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.

Ductal carcinoma in situ

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’ [380]. 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 [463]. 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 [232] 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 ([563] used the radially symmetric tumour model of Byrne and Chaplain [89] 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 [186188], 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 [175]). 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 [186]. 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 [187], 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 [186]. 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 [188]. 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 [357].

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, 224226, 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. [230]). An example of this is shown in figure 19 using a discrete cellular automaton model (described below).

Figure 19
Study of ductal carcinoma in situ (DCIS) using the mathematical model of Gatenby and co-workers [230] (colour online). Simulations show potential evolutionary pathways in carcinoma in situ. (a) Simulations start with a single layer of normal epithelial ...

3. Discrete modelling

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 [12], Moreira and Deutsch [381], Drasdo [159], Araujo and McElwain [34], Quaranta et al [439], Hatzikirou et al [260], Nagy [389], Abbott et al [1] Byrne et al [88], Fasano et al [182], Galle et al [213], Drasdo and Höhme [163], Thorne et al [519], Anderson et al [27], Deisboeck et al [150], Anderson and Quaranta [30], Quaranta et al [438] and Zhang et al [573]. 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 [412]. 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 [449] 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.

Figure 20
Discrete tumour modelling showing several (eukaryotic) cells from Rejniak [446]. The dots are cell–boundary points which connected by linear springs (thin lines) to model the elastic cell membrane. The interior circles denote cell nucleii. Cell–cell ...
Figure 21
The phases of cell proliferation from Rejniak [446]. (a) The cell begins the mitotic cycle (the interphase); (b) the anaphase—formation of daughter nucleii and an increase in cell volume; (c) the telophase—formation and pinchoff of contractile ...
Figure 22
Evolution of a 2D avascular tumour from Rejniak [446]. Grey: proliferating cells; white: quiescent cells and black: necrotic cells. Reprinted with permission from Rejniak 2007 J. Theor. Biol. 247 194. Copyright © 2007 Elsevier.

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 [157] and Anderson et al [31] 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 [129], Zheng et al [575], Li et al [339], Macklin and Lowengrub [350]) and cellular automaton models (e.g., Anderson [26], Anderson et al [32] 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 [33]). 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 [499] 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 [527] 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 [469], cell-differentiation [273] and cell-polarity [570]. Jiang et al [294] 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 [373375], by including a model of contact inhibition, and by Bauer et al [51] in the context of tumour-induced angiogenesis in a heterogeneous tumour microenviroment. Rubenstein and Kaufman [458] recently studied the role of ECM in glioma invasion using a version of the GGH model. Scianna et al [473] 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 [422] 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 [445], Li et al [339], Macklin and Lowengrub [350], Gerlee and Anderson [234, 235] and Anderson et al [31].

Another Monte Carlo based model that incorporates finite cell size was developed by Drasdo et al [164] 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 [376]. Note that in related work, Palsson and Othmer [403] modelled cells as deformable viscoelastic ellipsoids (see also Dallon and Othmer [144]).

Drasdo and Höhme [161] 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 [190] estimates were made for the biomechanical and kinetic parameters in the model. In later work, Drasdo and Höhme [162] 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 [215] extended the model to incorporate cell–substrate contact-dependent cell death (anoikis [500]) 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 [160] 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 [94] who presented further analysis of the continuum model. Radszuweit et al [440] 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 [266] investigated the effects of biomechanics and nutrient to determine their relative importance on the growth of cell populations. Ramis-Conde et al [444] 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 [443] extended this model to simulate tumour cell intravasation into blood vessels and investigated the role of cadherins in metastasis. Galle et al [214] 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 [277], Meineke et al [371] and Brodland and Veldhuis [77], 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 [304]. In Gevertz and Torquato [238], 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 [237] to simulate tumour growth in confined, heterogeneous environments.

Anderson et al [29] 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 [398] 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 [408] 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 [408]. 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 [486]) and Gatenby et al [230] 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 [472] 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 [360] 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 [158] developed a lattice-gas cellular automaton method for simulating the growth and size-saturation of avascular multicell tumour spheroids (see also the book [154] 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 [262] to study travelling fronts characterizing glioma cell invasion. Very recently, Hatzikirou et al [259] 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 [121] 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 [10] 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 [11]. Effects of nutrient heterogeneity as well as growth and invasion of cancerous tissue were investigated. This work was later extended by Owen et al [399] 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 [32] 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 [32], 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 [235] 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. [129]) 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 [236] 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 [415] 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.

Figure 23
Evolution of a 2D avascular tumour from the discrete model of Anderson [26]. The gradations of colour correspond to different cell phenotypes. The centre region contains necrotic cells. At late times, the outermost region contains the most aggressive ...

Bartha and Rieger [50] 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 [330] 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 [548] extended the model by incorporating the Fahraeus effect and refining the simulation of oxygen. Welter et al [548] 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 [549] recently extended the model to account for an arterio-venous network.

Mantzaris [361] 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 [451]. Mansury and Deisboeck [359] 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 [355] 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 [542] 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 [571] 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 [572] to incorporate mutations representing a simplified tumour progression pathway. Very recently, Quaranta et al [438] examined interactions between the tumour microenvironment and cancer cells at various scales using cell-based discrete mathematical models.

Ramis-Conde et al [442] 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 [245]. 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.

4. Hybrid modelling

4.1. Overview

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. [167]) and the equation-free method (see below).

4.2. Background

In an early study, Othmer and Stevens [398] 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 [476] 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 [177] 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 [176] 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 [309] (see also the review by Li et al [338]) 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 [70] proposed an equation-free approach to model the collective dynamics of populations of coupled, heterogeneous biological oscillators.

4.3. Coupled continuum–discrete descriptions of tumour cells

Recently, Kim et al [317] and Stolarska et al [498] 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 [144] 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.

Figure 24
Hybrid tumour modelling: evolution of a tumour spheroid in the absence of the outer gel from the study by Othmer and co-workers [317]. Necrosis is represented by the inner (white) region; it is enclosed by the continuum quiescent region and the outer ...

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 [55] and Frieboes et al (in preparation) [192] 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 [555].

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 [55]

Scd(x,t)=iSicd,whereSicd=CextractϕV(x,t)G(xxi),
(78)

where Cextract is a rate of extraction, ϕV is the volume fraction of viable tumour cells and G(xxi) is a function localized around the region xxi. 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 Vi(t)=tdetecttSicd(x,t)dΩdt. 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:

Sdc=iSidc,whereSidc(x,t)=Cdeposit(1ϕV)G(xxi),
(79)

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 Vr,i(t)=tdetecttSidc(x,t)dΩdt 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 [193] 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 [55] for human brain tumour specimens.

Figure 25
Evolution of a vascularized tumour using a hybrid continuum–discrete model for tumour cells. Discrete cells (blue) dots are released from hypoxic regions of the continuous tumour regions (grey). Discrete cells are converted back to continuum volume ...

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

E=12ijJij(x,y)ϕi(x)ϕj(y)dxdy,
(80)

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δ(xxj) 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

E=Ecc+Ecd+Edd,
(81)

where Ecc represents the self-interactions among the continuum fields

Ecc=12i,jNJij(x,y)ϕi(x)ϕj(y)dxdy,
(82)

the energy Edd represents the self-interactions among the discrete cells

Edd=12i,j>NJij(xi,xj)wiwj,
(83)

while the remaining energy describes the interactions between the discrete and continuous components

Ecd=12iNj>NwjJij(x,xj)ϕi(x)dx+12jNi<NwiJij(xi,y)ϕj(y)dy.
(84)

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:

δEδϕi=δEccδϕi+δEcdδϕi,
(85)

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(xy) is an even function and Jij = Jji, we obtain

δEccδϕi=jNJij(yx)ϕj(y)dy,
(86)
δEcdδϕi=j>NwjJij(xxj).
(87)

Assuming the kernel Jij is highly localized, a gradient approximation may be used for equation (86) to yield an energy of the form given in equation (53).

For the discrete cells, we may take the over-damped evolution

x.i=ζiδEδxi=ζi(δEddδxi+δEcdδxi),
(88)

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 [442] 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:

δEcdδxi=wijNJij(xxi)ϕj(x)dx.
(89)

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.

5. Conclusion

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 [57]. This is confirmed by increasing the cell–cell adhesion (or decreasing cell mobility) in mathematical models (see figure 12 and [126]).

Conversely, increasing nutrient levels in the tumour microenvironment may lead to greater morphological stability of tumours making it easier to remove surgically. This suggests [126] that treatments that seek to normalize tumour vasculature (by selectively ‘pruning’ weak blood vessels with targeted anti-angiogenic therapy [292]) 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 [405], Wicha et al [552], Dalerba et al [143], and Visvader and Lindeman [535]. 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 [172] and Ashkenazi et al [41]. Generic tumour models incorporating cancer stem cells have also been developed, e.g. see Ganguly and Puri [219], Piotrowska [416] and Galle et al [214]. Endering et al [173] 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 [174].

For further references on mathematical modelling of cancer stem cells we refer the reader to the recent review by Michor [378]. 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 [534] 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.

Acknowledgments

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.

References

1. Abbott RG, Forrest S, Pienta KJ. Simulating the hallmarks of cancer. Art. Life. 2006;12:617–34. [PubMed]
2. Abia LM, Angulo O, Lopez-Marcos JC. Age-structured population models and their numerical solution. Ecol. Modelling. 2005;188:112–36.
3. Adam JA. A simplified mathematical model of tumor growth. Math. Biosci. 1986;81:229–44.
4. Adam JA. A mathematical model of tumor growth: II. Effects of geometry and spatial nonuniformity on stability. Math. Biosci. 1987;86:183–211.
5. Adam JA. A mathematical model of tumor growth: III. Comparison with experiment. Math. Biosci. 1987;86:213–27.
6. Adam JA. General aspects of modeling tumor growth and the immune response. In: Adam J, Bellomo N, editors. A Survey of Models on Tumor Immune Systems Dynamics. Birkhauser; Boston, MA: 1996. pp. 15–87.
7. Adam JA, Bellomo NA. Survey of Models on Tumour Immune Systems Dynamics. Birkhauser; Boston, MA: 1997.
8. Adam JA, Maggelakis SA. Mathematical models of tumor growth: IV. Effects of a necrotic core. Math. Biosci. 1989;97:121–36. [PubMed]
9. Addison-Smith B, McElwain DLS, Maini PK. A simple mechanistic model of sprout spacing in tumour-associated angiogenesis. J. Theor. Biol. 2008;250:1–15. [PubMed]
10. Alarcón T, Byrne HM, Maini PK. A cellular automaton model for tumour growth in inhomogeneous environment. J. Theor. Biol. 2003;225:257–74. [PubMed]
11. Alarcón T, Byrne HM, Maini PK. A multiple scale model for tumor growth Multiscale. Model. Simul. 2005;3:440–75.
12. Alber MS, Kiskowski MA, Glazier JA, Jiang Y. In: On Cellular Automaton Approaches to Modeling Biological Cells (IMA Series on Mathematical Systems Theory in Biology, Communication and Finance. Rosenthal R, Gilliam DS, editors. Vol. 142. Springer; New York: 2002. pp. 1–40.
13. Ben Amar M, Goriely A. Growth and instability in soft tissues. J. Mech. Phys. Solids. 2005;53:2284–319.
14. Ambrosi D. Cellular traction as an inverse problem. SIAM J. Appl. Math. 2006;66:2049–60.
15. Ambrosi D, Bussolino F, Preziosi L. A review of vasculogenesis models. Comput. Math. Methods Med. 2005;6:1–19.
16. Ambrosi D, Duperray A, Peschetola V, Verdier C. Traction patterns of tumor cells. J. Math. Biol. 2009;58:163–81. [PubMed]
17. Ambrosi D, Gamba A, Serini G. Cell directional and chemotaxis in vascular morphogenesis. Bull. Math. Biol. 2004;66:1851–73. [PubMed]
18. Ambrosi D, Guana F. Mechanical aspects of growth in soft tissues. Boll. Unione Mat. Ital. 2004;7:775–81.
19. Ambrosi D, Guana F. Stress-modulated growth. Math. Mech. Solids. 2007;12:319–43.
20. Ambrosi D, Guillou ES, Di Martino ES. Stress modulated remodeling of a nonhomogeneous body. Biomech. Model. Mechanobiol. 2008;7:63–76. [PubMed]
21. Ambrosi D, Mollica F. On the mechanics of a growing tumor. Int. J. Eng. Sci. 2002;40:1297–16.
22. Ambrosi D, Mollica F. The role of stress in the growth of a multicell spheroid. J. Math. Biol. 2004;48:477–99. [PubMed]
23. Ambrosi D, Preziosi L. On the closure of mass balance models for tumor growth. Math. Models. Methods Appl. Sci. 2002;12:737–54.
24. Ambrosi D, Preziosi L. Cell adhesion mechanisms and elasto-viscoplastic mechanics of tumours. Mech. Model. Mechanobiol. 2009;8:397–413. [PubMed]
25. Andersen H, Mejlvang J, Mahmood S, Gromova I, Gromov P, Lukanidin E, Kriajevska M, Mellon JK, Tulchinsky E. Immediate and delayed effects of e-cadherin inhibition on gene regulation and cell motility in human epidermoid carcinoma cells. Mol. Cell. Biol. 2005;25:9138–50. [PMC free article] [PubMed]
26. Anderson ARA. A hybrid mathematical model of solid tumour invasion: the importance of cell adhesion. Math. Med. Biol. 2005;22:163–6. [PubMed]
27. Anderson ARA, Chaplain M, Rejniak K, Fozard J. Single-cell based models in biology and medicine. Math. Med. Biol. 2008;25:185–6.
28. Anderson ARA, Chaplain MAJ. Continuous and discrete mathematical models of tumor-induced angiogenesis. Bull. Math. Biol. 1998;60:857–900. [PubMed]
29. Anderson ARA, Chaplain MAJ, Newman EL, Steele RJC, Thompson AM. Mathematical modeling of tumour invasion and metastasis. J. Theor. Med. 2000;2:129–54.
30. Anderson ARA, Quaranta V. Integrative mathematical oncology. Nature Rev. Cancer. 2008;8:227–44. [PubMed]
31. Anderson ARA, Rejniak KA, Gerlee P, Quaranta V. Microenvironment driven invasion: a multiscale multimodel investigation. J. Math. Biol. 2009;58:579–624. [PubMed]
32. Anderson ARA, Weaver AM, Commmings PT, Quaranta V. Tumor morphology and phenotypic evolution driven by selective pressure from the microenvironment. Cell. 2006;127:905–15. [PubMed]
33. Anderson MP, Srolovitz DJ, Grest GS, Sahini PS. Computer simulation of grain growth: I. Kinetics. Acta Metall. 1984;32:783–91.
34. Araujo RP, McElwain DLS. A history of the study of solid tumour growth: the contribution of mathematical modelling. Bull. Math. Biol. 2004;66:1039–91. [PubMed]
35. Araujo RP, McElwain DLS. A linear-elastic model of anisotropic tumor growth. Eur. J. Appl. Math. 2004;15:365–84.
36. Araujo RP, McElwain DLS. New insights into vascular collapse and growth dynamics in solid tumors. J. Theor. Biol. 2004;228:335–46. [PubMed]
37. Araujo RP, McElwain DLS. A mixture theory for the genesis of residual stresses in growing tissues: I. A general formulation. SIAM J. Appl. Math. 2005;65:1261–84.
38. Araujo RP, McElwain DLS. A mixture theory for the genesis of residual stresses in growing tissues: II. Solutions to the biphasic equations for a multicell spheroid. SIAM J. Appl. Math. 2005;66:447–67.
39. Armstrong NJ, Painter KJ, Sherratt JA. A continuum approach to modeling cell–cell adhesion. J. Theor. Biol. 2006;243:98–113. [PMC free article] [PubMed]
40. Armstrong PB. Light and electron microscope studies of cell sorting in combinations of chick embryo and neural retina and retinal pigment epithelium. Wilhelm Roux Arch. 1971;168:125–41.
41. Ashkenazi R, Gentry SN, Jackson TL. Pathways to tumorigenesis-modeling mutation acquistion in stem cells and their progeny. Neoplasia. 2008;10:1170–82. [PMC free article] [PubMed]
42. Astanin S, Preziosi L. Multiphase Models of Tumour Growth. In: Bellomo N, et al., editors. (Selected Topics on Cancer Modelling: Genesis–Evolution–Immune Competition–Therapy) Birkhauser; Boston: 2007.
43. Astanin S, Tosin A. Mathematical model of tumour cord growth along the source of nutrient. Math. Model. Nat. Phenom. 2007;2:153–77.
44. Astanin S, Tosin L. Bi-population model of tumour growth and the warburg effect in tumour cords. J. Theor. Biol. 2009;258:578–90. [PubMed]
45. Augustin HG. Tubes, branches, and pillars: The many ways of forming a new vasculature. Circ. Res. 2001;89:645–47. [PubMed]
46. Ayati BP, Webb GF, Anderson ARA. Computational methods and results for structured multiscale models of tumor invasion. Multiscale Model. Simul. 2005;5:1–20.
47. Balding D, McElwain DLS. A mathematical model of tumor-induced capillary growth. J. Theor. Biol. 1985;114:53–73. [PubMed]
48. Bardelli A, Basile ML, Audero E, Giordano S, Wennström S, Ménard S, Comoglio PM, Ponzetto C. Concomitant activation of pathways downstream of grb2 and pi 3-kinase is required for met-mediated metastasis. Oncogene. 1999;18:1139–46. [PubMed]
49. Bartels U, Hawkins C, Jing M, Ho M, Dirks P, Rutka J, Stephens D, Bouffet E. Vascularity and angiogenesis as predictors of growth in optic pathway/hypothalamic gliomas. J. Neurosurg. 2006;104:314–20. [PubMed]
50. Bartha K, Rieger H. Vascular network remodeling via vessel cooption, regression and growth in tumors. J. Theor. Biol. 2006;241:903–18. [PubMed]
51. Bauer AL, Jackson TL, Jiang Y. A cell-based model exhibiting branching and anastomosis during tumor-induced angiogenesis. Biophys. J. 2007;92:3105–21. [PubMed]
52. Bauer TW, et al. Targeting of urokinase plasminogen activator receptor in human pancreatic carcinoma cells inhibits c-met- and insulin-like growth factor-1 receptor-mediated migration and invasion and orthotopic tumor growth in mice. Cancer Res. 2005;65:7775–81. [PubMed]
53. Bazaliy BV, Friedman A. A free boundary problem for an elliptic-parabolic system: Application to a model of tumor growth. Commun. Partial Diff. Eqns. 2003;28:517–60.
54. Bazaliy BV, Friedman A. Global existence and asymptotic stability for an ellipticparabolic free boundary problem: an application to a model of tumor growth. Indiana Math. J. 2003;52:1265–304.
55. Bearer EL, Lowengrub JS, Chuang YL, Frieboes HB, Jin F, Wise SM, Ferrari M, Agus DB, Cristini V. Multiparameter computational modeling of tumor invasion. Cancer Res. 2009;69:4493–501. [PMC free article] [PubMed]
56. Bell CD, Waizbard E. Variability of cell size in primary and metastatic human breast carcinoma. Invasion Metastasis. 1986;6:11–20. [PubMed]
57. Bello L, et al. Combinatorial administration of molecules that simultaneously inhibit angiogenesis and invasion leads to increased therapeutic efficacy in mouse models of malignant glioma. Clin. Cancer Res. 2004;10:4527–37. [PubMed]
58. Bellomo N, Bellouquid A, de Angelis E. The modelling of immune competition by generalised kinetic (Boltzmann) models: review and research perspectives. Math. Comput. Modelling. 2003;37:65–86.
59. Bellomo N, de Angelis E, Preziosi L. Multiscale modeling and mathematical problems related to tumor evolution and medical therapy. J. Theor. Med. 2003;5:111–36.
60. Bellomo N, Li NK, Maini PK. On the foundations of cancer modelling: selected topics, speculations, and perspectives. Math. Models Methods Appl. Sci. 2008;4:593–646.
61. Bellomo N, Preziosi L. Modelling and mathematical problems related to tumor evolution and its interaction with the immune system. Math. Comput. Modelling. 2000;32:413–542.
62. Bernsen HJJA, van der Kogel AJ. Antiangiogenic therapy in brain tumor models. J. Neuro-oncol. 1999;45:247–55. [PubMed]
63. Bertuzzi A, Fasano A, Gandolfi A. A free boundary problem with unilateral constraints describing the evolution of a tumor cord under the influence of cell killing agents. SIAM J. Math. Anal. 2005;36:882–915.
64. Bertuzzi A, Fasano A, Gandolfi A. A mathematical model for tumor cords incorporating the flow of an interstitial fluid. Math. Models Methods Appl. Sci. 2005;15:1735–77.
65. Bertuzzi A, Fasano A, Gandolfi A, Sinisgalli C. Atp production and necrosis formation in a tumour spheroid model. Math. Model. Natural Phenom. 2007;2:30–46.
66. Betteridge R, Owen MR, Byrne HM, Alarcón T, Maini PK. The impact of cell crowding and active cell movement on vascular tumour growth. Networks Heterogen. Media. 2006;1:515–35.
67. Blatz PJ, Ko WL. Application of finite elasticity theory to the deformation of rubbery materials. Trans. Soc. Rheol. 1962;6:223–51.
68. Bloemendal HJ, Logtenberg T, Voest EE. New strategies in anti-vascular cancer therapy. Eur. J. Clin. Invest. 1999;29:802–9. [PubMed]
69. Boccaccio C, Andò M, Tamagnone L, Bardelli A, Michieli P, Battistini C, Comoglio PM. Induction of epithelial tubules by growth factor hgf depends on the stat pathway. Nature. 1998;391:285–8. [PubMed]
70. Bold KA, Zou Y, Kevrekidis IG, Henson MA. An equation-free approach to analyzing heterogeneous cell population dynamics. J. Math. Biol. 2007;55:331–52. [PubMed]
71. Boushaba K, Levine HA, Nilsen-Hamilton M. A mathematical model for the regulation of tumor dormancy based on enzyme kinetics. Bull. Math. Biol. 2006;68:1495–526. [PubMed]
72. Brandt B, Junker R, Griwatz C, Heidl S, Brinkmann O, Semjonow A, Assmann G, Zänker KS. Isolation of prostate-derived single cells and cell clusters from human peripheral blood. Cancer Res. 1996;56:4556–61. [PubMed]
73. Breward CJW, Byrne HM, Lewis CE. The role of cell–cell interactions in a two-phase model for avascular tumour growth. J. Math. Biol. 2002;45:125–52. [PubMed]
74. Breward CJW, Byrne HM, Lewis CE. A multiphase model describing vascular tumor growth. Bull. Math. Biol. 2003;65:609–40. [PubMed]
75. Brizel DM, Scully SP, Harrelson JM, Layfield LJ, Bean JM, Prosnitz LR, Dewhirst MW. Tumor oxygenation predicts for the likelihood of distant metastases in human soft tissue sarcoma. Cancer Res. 1996;56:941–3. [PubMed]
76. Brizel DM, Sibley GS, Prosnitz LR, Scher RL, Dewhirst MW. Tumor hypoxia adversely affects the prognosis of carcinoma of the head and neck. Int. J. Radiat. Oncol. Biol. Phys. 1997;38:285–9. [PubMed]
77. Brodland GW, Veldhuis JH. Computer simulations of mitosis and interdependencies between orientation, cell shape and epithelia reshaping. J. Biomech. 2002;35:673–81. [PubMed]
78. Bueno H, Ercole G, Zumpano A. Asymptotic behaviour of quasi-stationary solutions of a nonlinear problem modelling the growth of tumours. Nonlinearity. 2005;18:1629–42.
79. Bullitt E, Zeng D, Gerig G, Aylward S, Joshi S, Smith JK, Lin W, Ewend MG. Vessel tortuosity and brain tumor malignance: a blinded study. Acad. Radiol. 2005;12:1232–40. [PMC free article] [PubMed]
80. Burgess PK, Kulesa PM, Murray JD, Alvord EC., Jr The interaction of growth rates and diffusion coefficients in a three dimensional mathematical model of gliomas. J. Neuropathol. Exp. Neurol. 1997;56:704–13. [PubMed]
81. Burton AC. Rate of growth of solid tumours as a problem of diffusion. Growth. 1966;30:157–76. [PubMed]
82. Bussolino F, Arese M, Audero E, Giraudo E, Marchiò S, Mitola S, Primo L, Serini G. Cancer Modelling and Simulation chapter 1 Biological Aspects of Tumour Angiogenesis. Chapman and Hall/CRC; London: 2003. pp. 1–22.
83. Byers SW, Sommers CL, Hoxter B, Mercurio AM, Tozeren A. Role of e-cadherin in the response of tumor cell aggregates to lymphatic, venous and arterial flow: measurement of cell–cell adhesion strength. J. Cell Sci. 1995;108:2053–64. [PubMed]
84. Byrne HM. The effect of time delays on the dynamics of avascular tumor growth. Math. Biosci. 1997;144:83–117. [PubMed]
85. Byrne HM. The importance of intercellular adhesion in the development of carcinomas. IMA J. Math. Med. Biol. 1997;14:305–23. [PubMed]
86. Byrne HM. A comparison of the roles of localized and nonlocalised growth factors in solid tumour growth. Math. Models Methods Appl. Sci. 1999;9:541–68.
87. Byrne HM. A weakly nonlinear analysis of a model of avascular solid tumour growth. J. Math. Biol. 1999;39:59–89. [PubMed]
88. Byrne HM, Alarcón T, Owen MR, Webb SW, Maini PK. Modeling aspects of cancer dynamics: A review. Phi. Trans. R. Soc. A. 2006;364:1563–78. [PubMed]
89. Byrne HM, Chaplain MAJ. Growth of nonnecrotic tumors in the presence and absence of inhibitors. Math. Biosci. 1995;130:151–81. [PubMed]
90. Byrne HM, Chaplain MAJ. Mathematical models for tumour angiogenesis: Numerical simulations and nonlinear wave solutions. Bull. Math. Biol. 1995;57:461–86. [PubMed]
91. Byrne HM, Chaplain MAJ. Growth of necrotic tumors in the presence and absence of inhibitors. Math. Biosci. 1996;135:187–216. [PubMed]
92. Byrne HM, Chaplain MAJ. Modelling the role of cell–cell adhesion in the growth and development of carcinomas. Math. Comput. Modelling. 1996;24:1–17.
93. Byrne HM, Chaplain MAJ. Free boundary value problems associated with the growth and development of multicellular spheroids. Eur. J. Appl. Math. 1997;8:639–58.
94. Byrne HM, Drasdo D. Individual-based and continuum models of growing cell populations: a comparison. J. Math. Biol. 2009;58:657–87. [PubMed]
95. Byrne HM, King JR, McElwain DLS, Preziosi L. A two-phase model of solid tumour growth. Appl. Math. Lett. 2003;16:567–73.
96. Byrne HM, Matthews P. Asymmetric growth of models of avascular solid tumors: exploiting symmetries. IMA J. Math. Appl. Med. Biol. 2002;19:1–29. [PubMed]
97. Byrne HM, Preziosi L. Modelling solid tumour growth using the theory of mixtures. Math. Med. Biol. 2003;20:341–66. [PubMed]
98. Cahn JW, Hilliard JE. Free energy of a nonuniform system: I. Interfacial free energy. J. Chem. Phys. 1958;28:258–67.
99. Cairns RA, Kalliomaki T, Hill RP. Acute (cyclic) hypoxia enhances spontaneous metastasis of kht murine tumors. Cancer Res. 2001;61:8903–08. [PubMed]
100. Capasso V, Morale D. Stochastic modelling of tumour-induced angiogenesis. J. Math. Biol. 2009;58:219–33. [PubMed]
101. Carmeliot P, Jain RK. Angiogenesis in cancer and other diseases. Nature. 2000;407:249–57. [PubMed]
102. Casciari JJ, Sotirchos SV, Sutherland RM. Variations in tumor cell growth rates and metabolism with oxygen concentration, glucose concentration, and extracellular pH. J. Cell. Physiol. 1992;151:386–94. [PubMed]
103. Castro M, Molina-Paris C, Deisboeck TS. Tumor growth instability and the onset of invasion. Phys. Rev. E. 2005;72:041907. [PubMed]
104. Chaplain MAJ. Reaction-diffusion prepatterning and its potential role in tumour invasion. J. Biol. Syst. 1995;3:929–36.
105. Chaplain MAJ. Avascular growth, angiogenesis and vascular growth in solid tumours: the mathematical modelling of the stages of tumour development. Mathl. Comput. Modelling. 1996;23:47–87.
106. Chaplain MAJ. Pattern formation in cancer. In: Chaplain MAJ, et al., editors. On Growth and Form: Spatio-Temporal Pattern Formation in Biology. Wiley; New York: 2000.
107. Chaplain MAJ, Benson DL, Maini PK. Nonlinear diffusion of a growth inhibitory factor in multicell spheroids. Math. Biosci. 1994;121:1–13. [PubMed]
108. Chaplain MAJ, Britton NF. On the concentration profile of a growth inhibitory factor in multicell spheroids. Math. Biosci. 1993;115:233–45. [PubMed]
109. Chaplain MAJ, Ganesh M, Graham IG. Spatio-temporal pattern formation on spherical surfaces: numerical simulation and application to solid tumour growth. J. Math. Biol. 2001;42:387–423. [PubMed]
110. Chaplain MAJ, Graziano L, Preziosi L. Mathematical modelling of the loss of tissue compression responsiveness and its role in solid tumor development. Math. Med. Biol. 2006;23:197–229. [PubMed]
111. Chaplain MAJ, Lolas G. Mathematical modeling of cancer cell invasion of tissue: the role of the urokinase plasminogen activation system. Math. Models Methods Appl. Sci. 2005;15:1685–734.
112. Chaplain MAJ, Lolas G. Mathematical modeling of cancer invasion of tissue: dynamic heterogeneity. Networks Heterogen. Media. 2006;1:399–439.
113. Chaplain MAJ, McDougall SR, Anderson ARA. Mathematical modeling of tumor-induced angiogenesis. Ann. Rev. Biomed. Eng. 2006;8:233–57. [PubMed]
114. Chaplain MAJ, Sleeman BD. Modelling the growth of solid tumours and incorporating a method for their classification using nonlinear elasticity theory. J. Math. Biol. 1993;31:431–79. [PubMed]
115. Chaplain MAJ, Stuart A. A model mechanism for the chemotactic response of endothelial cells to tumor angiogenesis factor. IMA J. Math. Appl. Med. Biol. 1993;10:149–68. [PubMed]
116. Chapman SJ, Shipley R, Jawad R. Multiscale modeling of fluid transport in tumors. Bull. Math. Biol. 2008;70:2334–57. [PubMed]
117. Chauviere A, Preziosi L. A mathematical framework to model migration of a cell population in the extracellular matrix. In: Chauviere A, et al., editors. Cell Mechanics: From Single Cell Scale-based Models to Multiscale Modeling. Chapman and Hall/CRC Press; Boca Raton, FL: 2009.
118. Chauviere A, Hillen T, Preziosi L. Modeling cell movement in anisotropic and heterogeneous tissues. Networks Heterogen. Media. 2007;2:333–57.
119. Chauviere A, Preziosi L, Hillen T. Modeling the motion of a cell population in the extracellular matrix. Discrete Contin. Dyn. Syst. B. 2007;(Suppl):250–59.
120. Chen CY, Byrne HM, King JR. The influence of growth-induced stress from the surrounding medium on the development of multicell spheroids. J. Math. Biol. 2001;43:191–220. [PubMed]
121. Chen WY, Annamreddy PR, Fan LT. Modeling growth of a heterogeneous tumor. J. Theor. Biol. 2003;221:205–27. [PubMed]
122. Chen X, Cui S, Friedman A. A hyperbolic free boundary problem modeling tumor growth: asymptotic behavior. Trans. Am. Math. Soc. 2005;357:4771–804.
123. Coleman BD, Noll W. Thermodynamics of elastic materials with conduction and viscosity. Arch. Ration. Mech. 1963;13:167–78.
124. Condeelis J, Singer RH, Segall JE. The great escape: when cancer cells hijack the genes for chemotaxis and motility. Annu. Rev. Cell Dev. Biol. 2005;21:695–718. [PubMed]
125. Coniglio A, deCandia A, DiTalia S, Gamba A. Percolation and burgers’ dynamics in a model of capillary formation. Phys. Rev. E. 2004;69:051910. [PubMed]
126. Cristini V, Frieboes HB, Gatenby R, Caserta S, Ferrari M, Sinek J. Morphologic instability and cancer invasion. Clin. Cancer Res. 2005;11:6772–9. [PubMed]
127. Cristini V, Frieboes HB, Li X, Lowengrub JS, Macklin P, Sanga S, Wise SM, Zheng X. Nonlinear modeling and simulation of tumor growth. In: Bellomo N, editor. Modelling and Simulation in Science, Engineering and Technology. Birkhauser; Boston: 2008.
128. Cristini V, Li X, Lowengrub JS, Wise SM. Nonlinear simulations of solid tumor growth using a mixture model: invasion and branching. J. Math. Biol. 2009;58:723–63. [PMC free article] [PubMed]
129. Cristini V, Lowengrub JS, Nie Q. Nonlinear simulation of tumor growth. J. Math. Biol. 2003;46:191–224. [PubMed]
130. Cruywagen GC, Woodward DE, Tracqui P, Bartoo GT, Murray JD, Alvord EC., Jr The modeling of diffusive tumors. J. Biol. Systems. 1995;3:937–945.
131. Cui S. Analysis of a mathematical model for the growth of tumors under the action of external inhibitors. J. Math. Biol. 2002;44:395–426. [PubMed]
132. Cui S. Analysis of a free boundary problem modeling tumor growth. Acta Math. Sin. 2005;21:1071–82.
133. Cui S. Lie group action and stability analysis of stationary solutions for a free boundary problem modelling tumor growth. J. diff. Equ. 2009;246:1845–82.
134. Cui S. Well-posedness of a multidimensional free boundary problem modelling the growth of nonnecrotic tumors. J. Funct. Anal. 2007;245:1–18.
135. Cui S. Asymptotic stability of the stationary solution for a hyperbolic free boundary problem modeling tumor growth. SIAM J. Math. Anal. 2008;40:1692–724.
136. Cui S, Escher J. Asymptotic behaviour of solutions of a multidimensional moving boundary problem modeling tumor growth. Commun. Partial Diff. Eqns. 2008;33:636–55.
137. Cui S, Escher J. Well-posedness and stability of a multi-dimensional tumor growth model. Arch. Ration. Mech. Anal. 2009;191:173–93.
138. Cui S, Friedman A. Analysis of a mathematical model of the effect of inhibitors on the growth of tumors. Math. Biosci. 2000;164:103–37. [PubMed]
139. Cui S, Friedman A. A free boundary problem for a singular system of differential equations: an application to a model of tumor growth. Trans. Am. Math. Soc. 2003;255:3537–90.
140. Cui S, Friedman A. Formation of necrotic cores in the growth of tumors: analytic results. Acta Mat. Sci. 2006;26:781–96.
141. Cui S, Wei X. Global existence for a parabolic-hyperbolic free boundary problem modelling tumor growth. Acta Math. Appl. Sin. 2005;21:597–614.
142. Cui S, Xu S. Analysis of mathematical models for the growth of tumors with time delays in cell proliferation. J. Math. Anal. Appl. 2007;336:523–41.
143. Dalerba P, Cho RW, Clarke MF. Cancer stem cells: Models and concepts. Ann. Rev. Med. 2007;58:267–84. [PubMed]
144. Dallon JC, Othmer HG. How cellular movement determines the collective force generated by the dictyostelium discoideum slug. J. Theor. Biol. 2004;231:299–306. [PubMed]
145. Date K, Matsumoto K, Kuba K, Shimura H, Tanaka M, Nakamura T. Inhibition of tumor growth and invasion by a four-kringle antagonist (HGF/NK4) for hepatocyte growth factor. Oncogene. 1998;17:3045–54. [PubMed]
146. de Angelis E, Preziosi L. Advection-diffusion models for solid tumour evolution in vivo and related free boundary problem. Math. Models Meth. Appl. Sci. 2000;10:379–407.
147. de Pillis LG, Radunskaya AE, Wiseman CL. A validated mathematical model of cell-mediated immune response to tumor growth. Cancer Res. 2005;65:7950–8. [PubMed]
148. Debnath J, Mills KR, Collins NL, Reginato MJ, Muthuswarmy SK, Brugge JS. The role of apoptosis in creating and maintaining luminal space within normal and oncogene-expressing mammary acini. Cell. 2002;111:29–40. [PubMed]
149. Deisboeck TS, Berens ME, Kansal AR, Torquato S, Stemmer-Rachamimov AO, Chiocca EA. Pattern of self-organization in tumour systems: complex growth dynamics in a novel brain tumour spherical model. Cell Proliferation. 2001;34:115–34. [PubMed]
150. Deisboeck TS, Zhang L, Yoon J, Costa J. In silico cancer modeling: is it ready for prime time? Nature Clin. Practice Oncol. 2009;6:34–42. [PMC free article] [PubMed]
151. DeJaeger K, Kavanagh MC, Hill R. Relationship of hypoxia to metastatic ability in rodent tumors. Br. J. Cancer. 2001;84:1280–5. [PMC free article] [PubMed]
152. Demuth T, Berens ME. Molecular mechanisms of glioma cell migration and invasion. J. Neuro-Oncol. 2004;70:217–28. [PubMed]
153. Derycke L, van Marck V, Depypere H, Bracke M. Molecular targets of growth, differentiation, tissue integrity and ectopic cell death in cancer cells. Cancer Biother. Radiopharm. 2005;20:579–88. [PubMed]
154. Deutsch A, Dormann S. Cellular Automaton Modeling of Biological Pattern Formation. Birkhäuser; New York: 2005.
155. Diaz JI, Tello JI. On the mathematical controllability in a simple growth tumors model by the internal localized action of inhibitors. Nonlinear Anal. 2003;4:109–25.
156. DiCarlo A, Quiligotti S. Growth and balance. Mech. Res. Commun. 2002;29:449–56.
157. Dillon R, Owen M, Painter K. A single-cell based model of multicellular growth using the immersed interface method. In: Khoo BC, et al., editors. Contemporary Mathematics: Moving Interface Problems and Applications in Fluid Dynamics. Vol. 466. AMS; Providence, RI: 2008. pp. 1–16.
158. Dormann S, Deutsch A. Modeling of self-organized avascular tumor growth with a hybrid cellular automaton. In Silico Biol. 2002;2:393–406. [PubMed]
159. Drasdo D. On selected individual-based approaches to the dynamics of multicellular systems. In: Alt W, et al., editors. Multiscale Modeling. Birkhauser; Basel: 2003.
160. Drasdo D. Coarse graining in simulated cell populations. Adv. Complex Syst. 2005;8:319–63.
161. Drasdo D, Höhme S. Individual-based approaches to birth and death in avascular tumors. Math. Comput. Modelling. 2003;37:1163–75.
162. Drasdo D, Höhme S. A single-scale-based model of tumor growth in vitro: monolayers and spheroids. Phys. Biol. 2005;2:133–47. [PubMed]
163. Drasdo D, Höhme S. On the role of physics in the growth and pattern of multicellular systems: What we learn from individual-cell based models? J. Stat. Phys. 2007;128:287–345.
164. Drasdo D, Kree R, McCaskill JS. Monte-carlo approach to tissue cell populations. Phys. Rev. E. 1995;52:6635–57. [PubMed]
165. Dyson J, Sanchez E, Villela-Bressan R, Webb G. An age and spatially structured model of tumor invasion and haptotaxis. Discrete Contin. Dyn. Syst. B. 2007;8:45–60.
166. Dyson J, Villella-Bressan R, Webb G. An age and spatially structured model of tumor invasion and haptotaxis ii. Math. Pop. Studies. 2008;15:73–95.
167. E W, Engquist B, Li XT, Ren WQ, Vanden-Eijnden E. Heterogeneous multiscale methods: A review. Commun. Math. Phys. 2007;2:367–450.
168. Eble JA, Haier J. Integrins in cancer treatment. Curr. Cancer Drug Targets. 2006;6:89–105. [PubMed]
169. Elliot CM, Luckhaus SA. generalized diffusion equation for phase separation of a multi-component mixture with interfacial free energy. Inst. Math. Appl. 1991 report 887.
170. Elvin P, Garner AP. Tumour invasion and metastasis: challenges facing drug discovery. Curr. Opin. Pharmacol. 2005;5:374–81. [PubMed]
171. Enam SA, Rosenblum ML, Edvardsen K. Role of extracellular matrix in tumor invasion: migration of glioma cells along fibronectinpositive mesenchymal cell processes. Neurosurgery. 1998;42:599–608. [PubMed]
172. Enderling H, Chaplain MA, Anderson ARA, Vaidya JS. A mathematical model of breast cancer development, local treatment and recurrence. J. Theor. Biol. 2007;246:245–59. [PubMed]
173. Enderling H, Hlatky L, Hahnfeldt P. Migration rules: tumours are conglomerates of self-metastases. Br. J. Cancer. 2009;100:1917–25. [PMC free article] [PubMed]
174. Enderling H, Park D, Hlatky L, Hahnfeldt P. The importance of spatial distribution of stemness and proliferation state in determining tumor radioresponse. Math. Model. Nat. Phenom. 2009;4:117–33.
175. Epstein M, Johnston CR. On the exact speed and amplitude of solitary waves in fluidfilled elastic tubes. Proc. R. Soc. A. 2001;457:1195–213.
176. Erban R, Kevrekidis IG, Othmer HG. An equation-free computational approach for extracting population-level behavior from individual-based models of biological dispersal. Physica D. 2006;215:1–24.
177. Erban R, Othmer HG. From individual to collective behavior in bacterial chemotaxis. SIAM J. Appl. Math. 2004;65:361–91.
178. Erler JT, Bennewith KL, Nicolau M, Dornhoefer N, Kong C, Le Q-T, Chi J-TA, Jeffrey SS, Giaccia AJ. Lysyl oxidase is essential for hypoxia-induced metastasis. Nature. 2006;440:1222–6. [PubMed]
179. Esteban MA, Maxwell PH. If, a missing link between metabolism and cancer. Nature Med. 2005;11:1047–8. [PubMed]
180. Eustace BK, et al. Functional proteomic screens reveal an essential extracellular role for hs90 alpha in cancer cell invasiveness. Nature Cell. Biol. 2004;6:507–14. [PubMed]
181. Fang JS, Gillies RD, Gatenby RA. Adaptation to hypoxia and acidosis in carcinogenesis and tumor progression. Semin. Cancer Biol. 2008;18:330–37. [PMC free article] [PubMed]
182. Fasano A, Bertuzzi A, Gandolfi A. Complex Systems in Biomedicine. Springer; Milan: 2006. Mathematical modelling of tumour growth and treatment; pp. 71–108.
183. Fedkiw RP, Aslam T, Merriman B, Osher S. A non-oscillatory eulerian approach to interfaces in multimaterial flows (the ghost fluid method) J. Comput. Phys. 1999;152:457–92.
184. Festjens N, Vanden Berghe T, Vandenabeele P. Necrosis, a well orchestrated form of cell demise: signaling cascades, important mediators and concomitant immune response. Biochim. Biophys. Acta. 2006;1757:1371–87. [PubMed]
185. Forsythe JA, Jiang BH, Iyer NV, Leung SW, Agani F, Koos RD, Semenza GL. Activation of vascular endothelial growth factor gene transcription by hypoxia-inducible factor 1. Mol. Cell. Biol. 1996;16:4604–13. [PMC free article] [PubMed]
186. Franks SJ, Byrne HM, King JR, Underwood JCE, Lewis CE. Modeling the early growth of ductal carcinoma in situ of the breast. J. Math. Biol. 2003;47:424–52. [PubMed]
187. Franks SJ, Byrne HM, Mudhar HS, Underwood JCE, Lewis CE. Mathematical modeling of comedo ductal carcinoma in situ of the breast. Math. Med. Biol. 2003;20:277–308. [PubMed]
188. Franks SJ, Byrne HM, Underwood JCE, Lewis CE. Biological inferences from a mathematical model of comedo duct carcinoma in situ of the breast. J. Theor. Biol. 2005;232:523–43. [PubMed]
189. Franks SJ, King JR. Interactions between a uniformly proliferating tumor and its surrounding uniform material properties. Math. Med. Biol. 2003;20:47–89. [PubMed]
190. Freyer JP, Sutherland RM. A reduction in the in situ rates of oxygen and glucose consumption of cells in emt6/ro spheroids during growth. J. Cell. Physiol. 1985;124:516–24. [PubMed]
191. Frieboes HB, Edgerton ME, Fruehauf JP, Rose FAJ, Worrall LK, Gatenby RA, Ferrari M, Cristini V. Prediction of drug response in breast cancer using integrative experimental/computational modeling. Cancer Res. 2009;69:4484–92. [PMC free article] [PubMed]
192. Frieboes HB, Fang J, Chuang Y-L, Wise SM, Lowengrub JS, Cristini V. Nonlinear simulations of three-dimensional multispecies tumor growth: II. Tumor invasion and angiogenesis. J. Theor. Biol. submitted. [PMC free article] [PubMed]
193. Frieboes HB, Lowengrub JS, Wise SM, Zheng X, Macklin P, Bearer EL, Cristini V. Computer simulation of glioma growth and morphology. NeuroImage. 2007;37:S59–70. [PMC free article] [PubMed]
194. Frieboes HB, Zheng X, Sun C-H, Tromberg BJ, Gatenby R, Cristini V. An integrated computational/experimental model of tumor invasion. Cancer Res. 2006;66:1597–604. [PubMed]
195. Friedl P. Prespecification and plasticity: shifting mechanisms of cell migration. Curr. Opin. Cell Biol. 2004;16:14–23. [PubMed]
196. Friedl P, Hegerfeldt Y, Tilisch M. Collective cell migration in morphogenesis and cancer. Int. J. Dev. Biol. 2004;48:441–9. [PubMed]
197. Friedl P, Noble PB, Walton PA, Laird DW, Chauvin PJ, Tabah RJ, Black M, Zänker KS. Migration of coordinated cell clusters in mesenchymal and epithelial cancer explants in vitro. Cancer Res. 1995;55:4557–60. [PubMed]
198. Friedl P, Wolf A. Tumor cell invasion and migration: diversity and escape mechanisms. Nature Rev. Cancer. 2003;3:362–74. [PubMed]
199. Friedman A. A hierarchy of cancer models and their mathematical challenges. Discrete Cont. Dyn. Syst. Ser. B. 2004;4:147–59.
200. Friedman A. A free boundary problem for a coupled system of elliptic, hyperbolic and stokes equations. Int. Free Bound. 2006;8:247–61.
201. Friedman A, Bellomo N, Maini PK. Mathematical analysis and challenges arising from models of tumor growth. Math. Models Methods Appl. Sci. 2007;17:1751–72.
202. Friedman A, Hu B. Asymptotic stability for a free boundary problem arising in a tumor model. J. Diff. Eqns. 2005;227:598–639.
203. Friedman A, Hu B. Bifurcation from stability to instability for a free boundary problem arising in a tumor model. Arch. Ration. Mech. Anal. 2006;180:293–330.
204. Friedman A, Hu B. Bifurcation from stability to instability for a free boundary problem modeling tumor growth by stokes equation. J. Math. Anal. Appl. 2007;327:643–64.
205. Friedman A, Hu B. Stability and instability of Liapunov-Schmidt and Hopf bifurcation for a free boundary problem arising in a tumor model. Trans. Am. Math. Soc. 2008;360:5291–342.
206. Friedman A, Reitich F. Analysis of a mathematical model for the growth of tumors. J. Math. Biol. 1999;38:262–84. [PubMed]
207. Friedman A, Reitich F. On the existence of spatially patterned dormant malignancies in a model for the growth of non-necrotic vascular tumors. Math. Models Methods Appl. Sci. 2001;11:601–25.
208. Friedman A, Reitich F. Symmetry-breaking bifurcation of analytic solutions to free boundary problems: an application to a model of tumor growth. Trans. Am. Math. Soc. 2001;353:1587–634.
209. Fukumura D, Jain RK. Tumor microenvironment abnormalities: causes, consequences, and strategies to normalize. J. Cell. Biochem. 2007;101:937–49. [PubMed]
210. Fung YC. Biomechanics: motion, flow, stress and growth. Springer; New York: 1990.
211. Fung YC. Biomechanics: Material Properties of Living Tissues. Springer; New York: 1993.
212. Galaris D, Barbouti A, Korantzopoulos P. Oxidative stress in hepatic ischemiareperfusion injury: the role of antioxidants and iron chelating compounds. Curr. Pharma. Des. 2006;12:2875–90. [PubMed]
213. Galle J, Aust G, Schaller G, Beyer T, Drasdo D. Individual cell-based models of the spatial temporal organization of multicellular systems-achievements and limitations. Cytometry. 2006;69 A:704–10. [PubMed]
214. Galle J, Hoffmann M, Aust G. From single cells to tissue architecture-a bottom-up approach to modeling the spatio-temporal organization of complex multicellular systems. J. Math. Biol. 2009;58:261–83. [PubMed]
215. Galle J, Loeffler M, Drasdo D. Modeling the effect of deregulated proliferation and apoptosis on the growth dynamics of epithelial cell populations in vitro. Biophys. J. 2005;88:62–75. [PubMed]
216. Galle J, Preziosi L, Tosin A. Contact inhibition of growth described by a multiphase and an individual cell-based model. Appl. Math. Lett. 2009;22:1483–90.
217. Gamba A, Ambrosi D, Coniglio A, deCandia A, DiTalia S, Giraudo E, Serini G, Preziosi L, Bussolino F. Percolation, morphogenesis and burgers dynamics in blood vessels formation. Phys. Rev. Lett. 2003;90:118101. [PubMed]
218. Ganghoffer JF. Some issues related to growth and goal functions for continuum biological systems. Phil. Mag. 2005;85:4353–91.
219. Ganguly R, Puri IK. Mathematical model for the cancer stem cell hypothesis. Cell Proliferation. 2006;39:3–14. [PubMed]
220. Garbey M, Zouridakis G. Modeling tumor growth: from differential deformable models to growth prediction of tumors detected in pet images. Eng. Med. Biol. Soc. 2003;3:2687–90.
221. Garikipati K, Arruda EM, Grosh K, Narayanan H, Calve S. A continuum treatment of growth in biological tissues. J. Mech. Phys. Solids. 2004;52:1595–625.
222. Garner L, Lau YY, Jackson TL, Uhler MD, Jordan DW, Gilgenbach RM. Incorporating spatial dependence into a multicellular tumor spheroid growth model. J. Appl. Phys. 2005;98:124701.
223. Gatenby RA, Frieden BR. Application of information theory and extreme physical information to carcinogenesis. Cancer Res. 2002;62:3675–2684. [PubMed]
224. Gatenby RA, Gawlinski ET. A reaction–diffusion model of cancer invasion. Cancer Res. 1996;56:5745–53. [PubMed]
225. Gatenby RA, Gawlinski ET. The glycolytic phenotype in carcinogenesis and tumor invasion: insights through mathematical models. Cancer Res. 2003;63:3847–54. [PubMed]
226. Gatenby RA, Gawlinski ET, Gmitro AF, Kaylor B, Gillies RJ. Acid-mediated tumor invasion: a multidisciplinary study. Cancer Res. 2006;66:5216–5223. [PubMed]
227. Gatenby RA, Gawlinsky ET. The Tumour Microenvironment. Wiley; London: 2003. Mathematical models of tumour invasion mediated by transformation-induced alteration of microenvironmental pH; pp. 85–99. [PubMed]
228. Gatenby RA, Gillies RJA. microenvironmental model of carcinogenesis. Nature Rev. Cancer. 2008;8:56–61. [PubMed]
229. Gatenby RA, Maini PK, Gawlinski ET. Analysis of tumor as an inverse problem provides a novel theoretical framework for understanding tumor biology and therapy. Appl. Math. Lett. 2002;6215:339–45.
230. Gatenby RA, Smallbone K, Maini PK, Rose F, Averill J, Nagle RB, Worrall L, Gillies RJ. Cellular adaptations to hypoxia and acidosis during somatic evolution of breast cancer. Br. J. Cancer. 2007;97:646–53. [PMC free article] [PubMed]
231. Gatenby RA, Vincent TL. An evolutionary model of carcinogenesis. Cancer Res. 2003;63:6212–20. [PubMed]
232. Gavaghan DJ, Brady JM, Behrenbruch CP, Highnam RP, Maini PK. Breast cancer: Modelling and detection. J. Theor. Medicine. 2002;4:3–20.
233. Gerisch A, Chaplain MA. Mathematical modelling of cancer cell invasion of tissue: local and non-local models and the effect of adhesion. J. Theor. Biol. 2008;250:684–704. [PubMed]
234. Gerlee P, Anderson ARA. An evolutionary hybrid cellular automaton model of solid tumor growth. J. Theor. Biol. 2007;246:583–603. [PMC free article] [PubMed]
235. Gerlee P, Anderson ARA. Stability analysis of a hybrid cellular automaton model of cell colony growth. Phys. Rev. E. 2007;75:051911. [PubMed]
236. Gerlee P, Anderson ARA. A hybrid cellular automaton model of clonal evolution in cancer: The emergence of the glycolytic phenotype. J. Theor. Biol. 2008;250:705–22. [PMC free article] [PubMed]
237. Gevertz JL, Gillies GT, Torquato S. Simulating tumor growth in confined heterogeneous environments. Phys. Biol. 2008;5:036010. [PubMed]
238. Gevertz JL, Torquato S. Modeling the effects of vasculature evolution on early brain tumor growth. J. Theor. Biol. 2006;243:517–31. [PubMed]
239. Gillies RJ, Gatenby RA. Adaptive landscapes and emergent phenotypes: why do cancers have high glycolysis? J. Bioenerg. Biomem. 2007;39:251–7. [PubMed]
240. Gillies RJ, Gatenby RA. Hypoxia and adaptive landscapes in the evolution of carcinogenesis. Cancer Metastasis Rev. 2007;26:311–17. [PubMed]
241. Gillies RJ, Robey I, Gatenby RA. Causes and consequences of increased glucose metabolism of cancers. J. Nucl. Med. 2008;49:24S–42S. [PubMed]
242. Gimbrone MA, Cotran RS, Leapman SB, Folkman J. Tumor growth and neovascularization: an experimental model using the rabbit cornea. J. Nat. Cancer Inst. 1974;52:413–27. [PubMed]
243. Glazier JA, Garner F. Simulation of the differential adhesion driven rearrangement of biological cells. Phys. Rev. E. 1993;47:2128–54. [PubMed]
244. Glimm J, Marchesin D, McBryan O. A numerical method for 2 phase flow with an unstable interface. J. Comput. Phys. 1981;39:179–200.
245. Gomez-Mourelo P, Sanchez E, Casasus L, Webb GF. A fully continuous individual-based model of tumor cell evolution. C. R. Biol. 2008;331:823–36. [PubMed]
246. Graeber TG, et al. Hypoxia-mediated selection of cells with diminished apoptotic potential in solid tumors. Nature. 1996;379:88–91. [PubMed]
247. Graner F, Glazier JA. Simulation of biological cell sorting using a two-dimensional extended potts model. Phys. Rev. Lett. 1992;69:2013–16. [PubMed]
248. Graziano L, Preziosi L. Mechanics in tumor growth. In: Mollica F, et al., editors. Modeling of Biological Materials. Birkhauser; New York: 2007. pp. 267–328.
249. Green JEF, Waters SL, Shakesheff KM, Byrne HM. A mathematical model of liver cell aggregation in vitro. Bull. Math. Biol. 2009;71:906–30. [PubMed]
250. Greenspan HP. Models for the growth of a solid tumor by diffusion. Stud. Appl. Math. 1972;51:317–40.
251. Greenspan HP. On the growth and stability of cell cultures and solid tumors. J. Theor. Biol. 1976;56:229–42. [PubMed]
252. Groebe K, Erz S, Mueller-Klieser W. Glucose diffusion coefficients determined from concentration profiles in emt6 tumor spheroids incubated in radioactively labeled l-glucose. Adv. Exp. Med. Biol. 1994;361:619–25. [PubMed]
253. Guiot C, Del Santo PP, Deisboeck TS. Morphological instability and cancer invasion: a ‘splashing water drop’ analogy. Theor. Biol. Med. Model. 2007;4:4. [PMC free article] [PubMed]
254. Habib S, Molina-Paris C, Deisboeck TS. Complex dynamics of tumors: modeling an emerging brain tumor system with coupled reaction-diffusion equations. Phys. A Stat. Mech. Appl. 2003;327:501–24.
255. Hanahan D, Weinberg RA. The hallmarks of cancer. Cell. 2000;100:57–70. [PubMed]
256. Harpold HL, Alvord EC, Jr, Swanson KR. The evolution of mathematical modeling of glioma proliferation and invasion. J. Neuropathol. Exp. Nerol. 2007;66:1–9. [PubMed]
257. Harris AL. Hypoxia-a key regulatory factor in tumor growth. Nature Rev. Cancer. 2002;2:38–47. [PubMed]
258. Hashizume H, Baluk P, Morikawa S, McLean JW, Thurston G, Roberge S, Jain RK, McDonald DM. Openings between defective endothelial cells explain tumor vessel leakiness. Am. J. Pathol. 2000;156:1363–80. [PubMed]
259. Hatzikirou H, Brusch L, Deutsch A. From cellular automaton rules to an effective macroscopic mean field description. Acta Phys. Polonica B. 2009 at press.
260. Hatzikirou H, Deutsch A, Schaller C, Simon M, Swanson K. Mathematical modeling of glioblastoma tumour development: A review. Math. Models Methods Appl. Sci. 2005;15:1779–94.
261. Hayot C, Debeir O, Van Ham P, Van Damme M, Kiss R, Decaestecker C. Characterization of the activities of actin-affecting drugs on tumor cell migration. Toxicol. Appl. Pharma. 2006;211:30–40. [PubMed]
262. Hatzikirou H, Brusch L, Schaller C, Simon M, Deustsch A. Prediction of traveling front behaviour in a lattice-gas cellular automaton model for tumour invasion. Comp. Appl. Math. 2009 at press.
263. Höckel M, Schlenger K, Aral B, Mitze M, Schäffer U, Vaupel P. Association between tumor hypoxia and malignant progression in advanced cancer of the uterine cervix. Cancer Res. 1996;56:4509–15. [PubMed]
264. Höckel M, Schlenger K, Höckel S, Vaupel P. Hypoxic cervical cancers with low apoptotic index are highly aggressive. Cancer Res. 1999;59:4525–8. [PubMed]
265. Höckel M, Schlenger K, Vaupel P. Tumor hypoxia: definitions and current clinical, biologic, and molecular aspects. J. Natl Cancer Inst. 2001;93:266–76. [PubMed]
266. Höhme S, Drasdo D. Biomechanical versus nutrient control: What determines the growth dynamics of mammalian cell populations? Math. Pop. Stud. 2009 at press.
267. Hegerfeldt Y, Tusch M, Bröcker EB, Friedl P. Collective cell movement in primary melanoma explants: plasticity of cell–cell interaction, β 1-integrin function, and migration strategies. Cancer Res. 2002;62:2125–30. [PubMed]
268. Helmlinger G, Netti PA, Lichtenbeld HC, Melder RJ, Jain RK. Solid stress inhibits the growth of multicellular tumor spheroids. Nat. Biotech. 1997;15:778–83. [PubMed]
269. Helmlinger G, Yuan F, Dellian M, Jain RK. Interstitial ph and po2 gradients in solid tumors in vivo: high-resolution measurements reveal a lack of correlation. Nature Med. 1997;3:177–82. [PubMed]
270. Hillen T, Painter K. A user's guide to pde models for chemotaxis. J. Math. Biol. 2009;58:183–217. [PubMed]
271. Hillen T, Painter K, Schmeiser C. Global existence for chemotaxis with finite sampling radius. Discrete Contin. Dyn. Syst. B. 2007;7:125–44.
272. Hogea CS, Murray BT, Sethian JA. Simulating complex tumor dynamics from avascular to vascular growth using a general level-set method. J. Math. Biol. 2006;53:86–134. [PubMed]
273. Hogeweg P. Evolving mechanisms of morphogenesis: on the interplay between differential adhesion and cell-differentiation. J. Theor. Biol. 2000;203:317–33. [PubMed]
274. Holash J, Maisonpierre PC, Compton D, Boland P, Alexander CR, Zagzag D, Yancopoulos GD, Wiegand SJ. Vessel cooption, regression, and growth in tumors mediated by angiopoietins and VEGF. Science. 1999;284:1994–8. [PubMed]
275. Holash J, Wiegand SJ, Yancopoulos GD. New model of tumor angiogenesis: dynamic balance between vessel regression and growth mediated by angiopoietins and VEGF. Oncogene. 1999;18:5356–62. [PubMed]
276. Holmes M, Sleeman B. A mathematical model of tumor angiogenesis incorporating cellular traction and viscoelastic effects. J. Theor. Biol. 2000;202:95–112. [PubMed]
277. Honda H, Tenemu M, Yoshida A. Differentiation of wing epidermal scale cells in a butterflu under the lateral inhibition model-appearance of large cells in a polygonal pattern. Acta Biotheor. 2000;48:121–36. [PubMed]
278. Huang Q, Shen HM, Ong CN. Emodin inhibits tumor cell migration through suppression of the phosphatidylinositol 3-kinase-cdc42/rac1 pathway. Cell. Mol. Life Sci. 2005;62:1167–75. [PubMed]
279. Humphrey JD. Continuum biomechanics of soft biological tissues. Proc. Roy. Soc. Lond A. 2003;459:303–11.
280. Humphrey JD, Rajagopal KR. A constrained mixture model for growth and remodeling of soft tissues. Math. Mod. Methods Appl. Sci. 2002;12:407–30.
281. Hurwitz H, et al. Bevacizumab plus irinotecan, fluorouacil, and leucovorin for metastatic colorectal cancer. New Eng. J. Med. 2004;350:2335–42. [PubMed]
282. Rohzin J, Sameni M, Ziegler G, Sloane BF. Pericellular ph affects distribution and secretion of cathepsin b in malignant cells. Cancer Res. 1994;54:6517–625. [PubMed]
283. Jackson TL. Intracellular accumulation and mechanism of action of doxorubicin in a spatio-temporal tumor model. J. Theor. Biol. 2003;220:201–13. [PubMed]
284. Jackson TL. A mathematical investigation of the multiple pathways to recurrent prostate cancer: Comparison with experimental data. Neoplasia. 2004;6:697–704. [PMC free article] [PubMed]
285. Jackson TL. A mathematical model of prostate tumor growth and androgenindependent relapse. Discrete Contin. Dyn. Syst. B. 2004;4:187–201.
286. Jackson TL, Ashkenazi R, Heusel S, Jain HV. Cancer modeling: A perspective on what's new and what's next. Contemp. Math. 2006;40:229–34.
287. Jackson TL, Byrne HM. A mechanical model of tumor encapsulation and transcapsular spread. Math. Biosci. 2002;180:307–28. [PubMed]
288. Jain HV, Nor JE, Jackson TL. Modeling the vegf-bcl-2-cxcl8 pathway in intratumoral angiogenesis. Bull. Math. Biol. 2008;70:89–117. [PubMed]
289. Jain RK. Determinants of tumor blood flow: a review. Cancer Res. 1988;48:2641–58. [PubMed]
290. Jain RK. Physiological barriers to delivery of monoclonal antibodies and other macromolecules in tumors. Cancer Res. 1990;50:814s–19s. [PubMed]
291. Jain RK. Delivery of molecular medicine to solid tumors: lessons from in vivo imaging of gene expression and function. J. Control. Release. 2001;74:7–25. [PubMed]
292. Jain RK. Normalizing tumor vasculature with anti-angiogenic therapy: a new paradigm for combination therapy. Nature Med. 2001;7:987–89. [PubMed]
293. Jbabdi S, Mandonnet E, Duffau H, Capelle L, Swanson KR, Pelegrini-Isaac M, Guillevin Remy, Benali H. Simulation of anisotropic growth of low-grade gliomas using diffusion tensor imaging. Mag. Res. Med. 2005;54:616–24. [PubMed]
294. Jiang Y, Pjesivac-Grbovic J, Cantrell C, Freyer JP. A multiscale model for avascular tumor growth. Biophys. J. 2005;89:3884–94. [PubMed]
295. Johnston MD, Edwards CM, Bodmer WF, Maini PK, Chapman SJ. Examples of mathematical modeling: tales from the crypt. Cell Cycle. 2007;6:2106–112. [PMC free article] [PubMed]
296. Johnston MD, Edwards CM, Bodmer WF, Maini PK, Chapman SJ. Mathematical modeling of cell population dynamics in the colonic crypt and in colorectal cancer. Proc. Natl. Acad. USA. 104:4008–13. [PubMed]
297. Jolly C, Morimoto RI. Role of the heat shock response and molecular chaperones in oncogenesis and cell death. J. Natl. Cancer Inst. 2000;92:1564–72. [PubMed]
298. Jones AF, Byrne HM, Gibson JS, Dold JW. Mathematical model for the stress induced during avascular tumor growth. J. Math. Biol. 2000;40:473–99. [PubMed]
299. Jones PF, Sleeman BD. Angiogenesis-understanding the mathematical challenge. Angiogenesis. 2006;9:127–38. [PubMed]
300. Joseph DD. Fluid Dynamics of Viscoelastic Liquids. Springer; New York: 1990.
301. Kallinowski F, Vaupel P, Runkel S, Berg G, Fortmeyer HP, Baessler KH, Wagner K, Mueller-Klieser W, Walenta S. Glucose uptake, lactate release, ketone body turnover, metabolic milieu and ph distributions in human cancer xenografts in nude rats. Cancer Res. 1988;48:7264–72. [PubMed]
302. Kaneko K, Satoh K, Masamune A. T. myosin light chain kinase inhibitors can block invasion and adhesion of human pancreatic cancer cell lines. Pancreas. 2002;24:34–41. [PubMed]
303. Kansal AA, Torquato S, Harsh GR, IV, Chiocca EA, Deisboeck TS. Cellular automaton of idealized brain tumor growth dynamics. BioSystems. 2000;55:119–27. [PubMed]
304. Kansal AA, Torquato S, Harsh GR, IV, Chiocca EA, Deisboeck TS. Emergence of a subpopulation in a computational model of tumor growth. J. Theor. Biol. 2000;207:431–41. [PubMed]
305. Kansal AA, Torquato S, Harsh GR, IV, Chiocca EA, Deisboeck TS. Simulated brain tumor growth dynamics using a 3-d cellular automaton. J. Theor. Biol. 2000;203:367–82. [PubMed]
306. Kaur B, Khwaja FW, Severson EA, Matheny SL, Brat DJ, VanMeir EG. Hypoxia and the hypoxia-inducible-factor pathway in glioma growth and angiogenesis. Neuro-Oncol. 2005;7:134–53. [PMC free article] [PubMed]
307. Keller PJ, Pampaloni F, Stelzer EHK. Life sciences require the third dimension. Curr. Opin. Cell Biol. 2006;18:117–24. [PubMed]
308. Kenny PA, Lee GY, Bissell MJ. Targeting the tumor microenvironment. Front. Biosci. 2007;12:3468–74. [PMC free article] [PubMed]
309. Kevrekidis IG, Gear CW, Hyman JM, Kevrekidis P, Runborg O, Theodoropoulos K. Equation-free, coarse-grained multiscale computation: enabling microscopic simulators to perform system-level analysis. Comm. Math. Sci. 2003;1:715–62.
310. Kevrekidis PG, Whitaker N, Good DJ, Herring GJ. Minimal model for tumor angiogenesis. Phys. Rev. E. 2006;73:061926. [PubMed]
311. Khain E, Sander LM. Dynamics of pattern formation in invasive tumor growth. Phys. Rev. Lett. 2006;96:188103. [PubMed]
312. Khain E, Sander LM. Generalized cahn-hilliard equation for biological applications. Phys. Rev. E. 2008;77:051129. [PubMed]
313. Khain E, Sander LM, Schneider-Mizell CM. The role of cell–cell adhesion in wound healing. J. Stat. Phys. 2007;128:209–18.
314. Khain E, Sander LM, Stein AM. A model for glioma growth. Complexity. 2005;11:53–7.
315. Kim JB. Three-dimensional tissue culture models in cancer biology. J. Biomol. Screening. 2005;15:365–77. [PubMed]
316. Kim JS, Lowengrub JS. Phase field modeling and simulation of three phase flows. Int. Free Bound. 2005;7:435–66.
317. Kim Y, Stolarska MA, Othmer HG. A hybrid model for tumor spheroid growth in vitro: I. Theoretical development and early results. Math. Methods Appl. Sci. 2007;17:1773–98.
318. Kloner RA, Jennings RB. Consequences of brief ischemia: stunning, preconditioning, and their clinical implications: Part I. Circulation. 2001;104:2981–9. [PubMed]
319. Kopfstein L, Christofori G. Metastasis: cell-autonomous mechanisms versus contributions by the tumor microenvironment. Cell. Mol. Life Sci. 2006;63:449–68. [PubMed]
320. Kuhl E, Holzapfel GA. A continuum model for remodeling in living structures. J. Mater. Sci. 2007;42:8811–23.
321. Kuhl E, Menzel A, Steinmann P. Computational modeling of growth: A critical review, a classification of concepts and two new consistent approaches. Comput. Mech. 2003;32:71–88.
322. Kuiper RAJ, Schellens JHM, Blijham GH, Beijnen JH, Voest EE. Clinical research on antiangiogenic therapy. Pharmacol. Res. 1998;37:1–16. [PubMed]
323. Kunkel P, Ulbricht U, Bohlen P, Fillbrandt R, Brockmann MA, Stavrou D, Westphal M, Lamszus K. Inhibition of glioma angiogenesis and growth in vivo by systemic treatment with a monoclonal antibody against vascular endothelial growth factor receptor-2. Cancer Res. 2001;61:6624–8. [PubMed]
324. Kunz-Schughart LA, Freyer JP, Hofstaedter F, Ebner R. The use of 3-d cultures for high-throughput screening: the multicellular spheroid model. J. Biomol. Screening. 2004;9:273–85. [PubMed]
325. Kuusela E, Alt W. Continuum model of cell adhesion and migration. J. Math. Biol. 2009;58:135–61. [PubMed]
326. Lah TT, Alonso MBD, van Noorden CJF. Antiprotease therapy in cancer: hot or not? Exp. Opin. Biol. Ther. 2006;6:257–79. [PubMed]
327. Lamszus K, Kunkel P, Westphal M. Invasion as limitation to anti-angiogenic glioma therapy. Acta Neurochir Suppl. 2003;88:169–77. [PubMed]
328. Landman KA, Please CP. Tumour dynamics and necrosis: Surface tension and stability. IMA J. Math. Appl. Med. Biol. 2001;18:131–58. [PubMed]
329. Lanza V, Ambrosi D, Preziosi L. Exogenous control of vascular network formation in vitro: a mathematical model. Networks Heterogen. Media. 2006;1:621–37.
330. Lee DS, Rieger H, Bartha K. Flow correlated percolation during vascular remodeling in growing tumors. Phys. Rev. Lett. 2006;96:058104. [PubMed]
331. Lester RD, Jo M, Campana WM, Gonias SL. Erythropoietin promotes mcf-7 breast cancer cell migration by an erk/mitogenactivated protein kinase-dependent pathway and is primarily responsible for the increase in migration observed in hypoxia. J. Biol. Chem. 2005;280:39273–7. [PubMed]
332. Levine H, Sleeman B. In: Modelling Tumour-Induced Angiogenesis Cancer Modelling and Simulation. Preziosi L, editor. Chapman&Hall/CRC; Boca Raton, FL: 2003. pp. 147–84.
333. Levine HA, Nilsen-Hamilton M. Angiogenesis-a biochemical/mathematical perspective. Tutor. Math. Biosci. III. 2006;1872:23–76.
334. Levine HA, Pamuk S, Sleeman BD, Nilsen-Hamilton M. Mathematical modeling of capillary formation and development in tumor angiogenesis: penetration into the stroma. Bull. Math. Biol. 2001;63:801–63. [PubMed]
335. Levine HA, Sleeman BD, Nilsen-Hamilton M. Mathematical modeling of the onset of capillary formation initiating angiogenesis. J. Math. Biol. 2001;42:195–238. [PubMed]
336. Levine HA, Smiley MW, Tucker AL, Nilsen-Hamilton MA. mathematical model for the onset of avascular tumor growth in response to the loss of p53 function. Cancer Informatics. 2006;2:163–88. [PMC free article] [PubMed]
337. Levine HA, Tucker AL, Nilsen-Hamilton MA. mathematical model for the role of cell signal transduction in the initiation and inhibition of angiogenesis. Growth Factors. 2002;20:155–75. [PubMed]
338. Li J, Kevrekidis PG, Gear CW, Kevrekidis IG. Deciding the nature of the coarse equation through microscopic simulations: The baby-bathwater scheme. SIAM Rev. 2007;49:469–87.
339. Li X, Cristini V, Nie Q, Lowengrub JS. Nonlinear three-dimensional simulation of solid tumor growth. Discrete Dyn. Contin. Dyn. Syst. B. 2007;7:581–604.
340. Lin IE, Tabor LA. A model for stress-induced growth in the developing heart. J. Biomech. Eng. 1995;117:343–9. [PubMed]
341. Liotta LA, Kohn EC. The microenvironment of the tumour-host interface. Nature. 2001;411:375–9. [PubMed]
342. Lloyd BA, Szczerba D, Rudin M, Szekely G. A computational framework for modeling solid tumour growth. Phil. Trans. R. Soc. A. 2008;366:3301–18. [PubMed]
343. Lloyd BA, Szczerba D, Szekely G. A coupled finite element model of tumor growth and vascularization. In: Ayache N, et al., editors. Medical Image Computing and Computer-Assisted Intervention-MICCA 2007 10th Int. Conf. Vol. 4792. Springer; New York: 2007. pp. 874–81. Lecture Notes in Computer Science. [PubMed]
344. Lockett J, Yin SP, Li XH, Meng YH, Sheng SJ. Tumor suppressive maspin and epithelial homeostasis. J. Cell. Biochem. 2006;97:651–60. [PubMed]
345. Lubkin SR, Jackson TL. Multiphase mechanics of capsule formation in tumors. J. Biomed. Eng.—Trans. ASME. 2002;124:237–43. [PubMed]
346. MacArthur BD, Please CP. Residual stress generation and necrosis formation in multi-cell tumour spheroids. J. Math. Biol. 2004;49:537–52. [PubMed]
347. Macklin P. Numerical Simulation of Tumor growth and Chemotherapy. University of Minnesota; Minnesota: 2003.
348. Macklin P, Lowengrub JS. Evolving interfaces via gradients of geometry-dependent interior Poisson problems: Application to tumor growth. J. Comput. Phys. 2005;203:191–220.
349. Macklin P, Lowengrub JS. An improved geometry-aware curvature discretization for level set methods: Application to tumor trowth. J. Comput. Phys. 2006;215:392–401.
350. Macklin P, Lowengrub JS. Nonlinear simulation of the effect of microenvironment on tumor growth. J. Theor. Biol. 2007;245:677–704. [PubMed]
351. Macklin P, Lowengrub JS. A new ghost cell/level set method for moving boundary problems: Application to tumor growth. J. Sci. Comput. 2008;35:266–99. [PMC free article] [PubMed]
352. Macklin P, McDougall S, Anderson ARA, Chaplain MAJ, Cristini V, Lowengrub J. Multiscale modeling and nonlinear simulation of vascular tumour growth. J. Math. Biol. 2009;58:765–98. [PMC free article] [PubMed]
353. Maggelakis SA, Adam JA. Mathematical model of prevascular growth of a spherical carcinoma. Math. Comput. Modelling. 1990;13:23–38.
354. Mallet DG, de Pillis LG. A cellular automata model of tumor-immune system interactions. J. Theor. Biol. 2006;239:334–50. [PubMed]
355. Mallett DG, Pettet GJ. A mathematical model of integrin-mediated haptotactic cell migration. Bull. Math. Biol. 2006;68:231–53. [PubMed]
356. Malvern LE. Introduction of the Mechanics of a Continuous Medium. Prentice Hall; NJ, Englewood Cliffs: 1969.
357. Man YG, Sang QXA. The significance of focal myoepithelial cell layer disruptions in human breast tumor invasion: a paradigm shift from the ‘protease-centered’ hypothesis. Exp. Cell Res. 2004;2:103–18. [PubMed]
358. Manoussaki D, Lubkin SR, Vernon RB, Murray JD. A mechanical model for the formation of vascular networks in vitro. Acta Biotheor. 1996;44:271–82. [PubMed]
359. Mansury Y, Deisboeck TS. Modeling tumors as complex biosystems: An agentbased approach. In: Deisboeck TS, Kresh JY, editors. Complex Systems Science in Biomedicine. Springer; New York: 2006. pp. 573–602.
360. Mansury Y, Kimura M, Lobo J, Deisboeck TS. Emerging patterns in tumor systems: simulating the dynamics of multicellular clusters with an agent-based spatial agglomeration model. J. Theor. Biol. 2002;219:343–70. [PubMed]
361. Mantzaris NV. Stochastic and deterministic simulations of heterogeneous cell population dynamics. J. Theor. Biol. 2006;241:690–706. [PubMed]
362. Mantzaris NV, Webb S, Othmer HG. Mathematical modeling of tumor-induced angiogenesis. J. Math. Biol. 2004;49:111–87. [PubMed]
363. Marchant BP, Norbury J, Byrne HM. Biphasic behaviour in malignant invasion. Math. Med. Biol. 2006;23:173–96. [PubMed]
364. Marchant BP, Norbury J, Sherratt JA. Travelling wave solutions to a haptotaxisdominated model of malignant invasion. Nonlinearity. 2001;14:1653–71.
365. Martinsa ML, Ferreira SC, Jr, Vilela MJ. Multiscale models for the growth of avascular tumors. Phys. Life Rev. 2007;4:128–56.
366. McDougall SR, Anderson ARA, Chaplain MAJ. Mathematical modeling of dynamic adaptive tumour-induced angiogenesis: Clinical applications and therapeutic targeting strategies. J. Theor. Biol. 2006;241:564–89. [PubMed]
367. McDougall SR, Anderson ARA, Chaplain MAJ, Sherratt J. Mathematical modelling of flow through vascular networks: implications for tumour-induced angiogenesis and chemotherapy strategies. Bull. Math. Biol. 2002;64:673–702. [PubMed]
368. McElwain DLS, Morris LE. Apoptosis as a volume loss mechanism in mathematical models of solid tumor growth. Math. Biosci. 1978;39:147–57.
369. McElwain DLS, Pettet GJ. Cell migration in multicell spheroids: swimming against the tide. Bull. Math. Biol. 1993;55:655–74. [PubMed]
370. McLean GW, Carragher NO, Avizienyte E, Evans J, Brunton VG, Frame MC. The role of focal-adhesion kinase in cancer: a new therapeutic opportunity. Nature Rev. Cancer. 2005;5:505–15. [PubMed]
371. Meineke FA, Potten CS, Loeffler M. Cell migration and organization in the intestinal crypt using a lattice-free model. Cell Proliferation. 2001;34:253–66. [PubMed]
372. Menzel A. Modelling of anisotropic growth in biological tissues- a new approach and computational aspects. Biomech. Model. Mechanobiol. 2005;3:147–71. [PubMed]
373. Merks RMH, Brodsky SV, Goligorksy MS, Newman SA, Glazier JA. Cell elongation is key to in silico replication of in vitro vasculogenesis and subsequent remodeling. Dev. Biol. 2006;289:44–54. [PMC free article] [PubMed]
374. Merks RMH, Glazier JA. Dynamic mechanisms of blood vessel growth. Nonlinearity. 2006;19:C1–10. [PMC free article] [PubMed]
375. Merks RMH, Perrynand ED, Shirinifard A, Glazier JA. Contact-inhibited chemotaxis in de novo and sprouting blood-vessel growth. PloS Comput. Biol. 2008;4:e1000163. [PMC free article] [PubMed]
376. Metropolis N, Rosenbluth AW, Rosenbluth MN, Teller AH, Teller E. Equation of state calculations by fast computing machines. J. Chem. Phys. 1953;21:1087–92.
377. Michieli P, Basilico C, Pennacchietti S, Maffè A, Tamagnone L, Giordano S, Bardelli A, Comoglio PM. Mutant met mediated transformation is liganddependent and can be inhibited by HGF antagonists. Oncogene. 1999;18:5221–31. [PubMed]
378. Michor F. Mathematical models of cancer stem cells. J. Clin. Oncol. 2008;26:2854–61. [PubMed]
379. Milde F, Bergdorf M, Koumoutsakos P. A hybrid model for three-dimensional simulations of sprouting angiogenesis. Biophys. J. 2008;95:3146–60. [PubMed]
380. Mills RR, Handy AM, Oberman HA. The breast. In: Stenberg SS, editor. Diagnostic Surgical Pathology. Vol. 1. Academic; New York: 1999. pp. 319–85.
381. Moreira J, Deutsch A. Cellular automaton models of tumor development: a critical review. Adv. Complex Syst. 2002;5:247–67.
382. Morotti A, Mila S, Accornero P, Tagliabue E, Ponzetto C. K252a inhibits the oncogenic properties of met, the HGF receptor. Oncogene. 2002;21:4885–93. [PubMed]
383. Mueller-Klieser W. Multicellular spheroids: a review on cellular aggregates in cancer research. J. Cancer Res. Clin. Oncol. 1987;113:101–22. [PubMed]
384. Mueller-Klieser W, Freyer JP, Sutherland RM. Influence of glucose and oxygen supply conditions on the oxygenation of multicellular spheroids. Br. J. Cancer. 1986;53:345–53. [PMC free article] [PubMed]
385. Mullins WW, Sekerka RF. Morphological instability of a particle growing by diffusion or heat flow. J. Appl. Phys. 1963;34:323–9.
386. Murray J, Oster G. Cell traction models for generation of pattern and form in morphogenesis. J. Math. Biol. 1984;33:489–520. [PubMed]
387. Murray JD. Interdisciplinary Applied Mathematics. 3rd edn Vol. 17. Springer; Berlin: 2002. Mathematical biology II: spatial models and biomedical applications.
388. Naganuma H, Kimurat R, Sasaki A, Fukamachi A, Nukui H, Tasaka K. Complete remission of recurrent glioblastoma multiforme following local infusions of lymphokine activated killer cells. Acta Neurochir. 1989;99:157–60. [PubMed]
389. Nagy JD. The ecology and evolutionary biology of cancer: A review of mathematical models of necrosis and tumor cell diversity. Math. Biosci. Eng. 2005;2:381–418. [PubMed]
390. Nelson CM, Bissell MJ. Of extracellular matrix, scaffolds, and signaling: tissue architecture regulates development, homeostasis, and cancer. Annu. Rev. Cell Dev. Biol. 2006;22:287–309. [PMC free article] [PubMed]
391. Neville AA, Matthews PC, Byrne HM. Interactions between pattern formation and domain growth. Bull. Math. Biol. 2006;68:1975–2003. [PubMed]
392. Ngwa G, Maini P. Spatio-temporal patterns in a mechanical model for mesenchymal morphogenesis. J. Math. Biol. 1995;33:489–520. [PubMed]
393. Nichols MG, Foster TH. Oxygen diffusion and reaction kinetics in the photodynamic therapy of multicell tumour spheroids. Phys. Med. Biol. 1994;39:2161–81. [PubMed]
394. Nör JE, Christensen J, Liu J, Peters M, Mooney DJ, Strieter RM, Polverini PJ. Up-regulation of BCL-2 in microvascular endothelial cells enhances intratumoral angiogenesis and accelerates tumor growth. Cancer Res. 2001;61:2183–8. [PubMed]
395. O'Connor JPB, Jackson A, Parker GJM, Jayson GC. Dce-mri biomarkers in the clinical evaluation of antiangiogenic and vascular disrupting agents. Br. J. Cancer. 2007;96:189–95. [PMC free article] [PubMed]
396. Orme ME, Chaplain MAJ. Two-dimensional models of tumour angiogenesis and anti-angiogenesis strategies. Math. Med. Biol. 1997;14:189–205. [PubMed]
397. Osher S, Sethian J. Fronts propagating with curvature-dependent speed: algorithms based on Hamilton–Jacobi formulation. J. Comput. Phys. 1988;79:12.
398. Othmer HG, Stevens A. Aggregration, blowup, and collapse: the abc's of taxis in reinforced random walks. SIAM. J. Appl. Math. 1997;57:1044–81.
399. Owen MR, Alarcón T, Maini PK, Byrne HM. Angiogenesis and vascular remodeling in normal and cancerous tissues. J. Math. Biol. 2009;58:689–721. [PubMed]
400. Padera TP, Stoll BR, Tooredman JB, Capen D, di Tomaso E, Jain R. Cancer cells compress intratumour vessels. Nature. 2004;427:695. [PubMed]
401. Page DL, Anderson TJ, Sakamoto G. Diagnostic Histopathology of the Breast. Churchill Livingstone; New York: 1987.
402. Paku S. First step of tumor-related angiogenesis. Lab. Invest. 1991;65:334–46. [PubMed]
403. Palsson E, Othmer HG. A model for individual and collective cell movement in dictyostelium discoideum. Proc. Natl Acad. Sci. USA. 2000;97:10338–453. [PubMed]
404. Pamuk S. Qualitative analysis of a mathematical model for capillary formation in tumor angiogenesis. Math. Models Methods Appl. Sci. 2003;13:19–33.
405. Pardal R, Clarke MF, Morrison SJ. Applying the principles of stem cell biology to cancer. Nature Rev. Cancer. 2003;3:895–902. [PubMed]
406. Parnuk S. A mathematical model for capillary formation and development in tumor angiogenesis: a review. Chemotherapy. 2006;52:35–7. [PubMed]
407. Patan S, Tanda S, Roberge S, Jones RC, Jain RK, Munn LL. Vascular morphogenesis and remodeling in a human tumor xenograft: blood vessel formation and growth after ovariectomy and tumor implantation. Circ. Res. 2001;89:732–9. [PMC free article] [PubMed]
408. Patel AA, Gawlinski ET, Lemieux SK, Gatenby RA. A cellular automaton model of early tumor growth and invasion: The effects of native tissue vascularity and increased anaerobic tumor metabolism. J. Theor. Biol. 2001;213:315–31. [PubMed]
409. Paweletz N, Knierim M. Tumor-related angiogenesis. Crit. Rev. Oncol. Hematol. 1989;9:197–242. [PubMed]
410. Peirce SM. Computational and mathematical modeling of angiogenesis. Microcirculation. 2008;15:739–51. [PMC free article] [PubMed]
411. Pennacchietti S, Michieli P, Galluzzo M, Giordano S, Comoglio P. Hypoxia promotes invasive growth by transcriptional activation of the met protooncogene. Cancer Cell. 2003;3:347–61. [PubMed]
412. Peskin C. The immersed boundary method. Acta Numer. 2002;11:479–517.
413. Peterson JW, Carey GF, Knezevic DJ, Murray BT. Adaptive finite element methodology for tumour angiogenesis modelling. Intl. J. Numer. Methods Eng. 2007;69:1212–38.
414. Pettet G, Please CP, Tindall M, McElwain D. The migration of cells in multicell tumor spheroids. Bull. Math. Biol. 2001;63:231–57. [PubMed]
415. Piotrowska MJ, Angus SD. A quantitative cellular automaton model of in vitro multicellular spheroid tumour growth. J. Theor. Biol. 2009;258:165–78. [PubMed]
416. Piotrowska MJ, Enderling H, van der Heiden U, Mackey MC. Mathematical modeling of stem cells related to cancer. In: Dittmar T, et al., editors. Complex Systems in Biomedicine. Nova; New York: 2008.
417. Plank MJ, Sleeman BD. A reinforced random walk model of tumour angiogenesis and anti-angiogenic strategies. Math. Med. Biol. 2003;20:135–81. [PubMed]
418. Plank MJ, Sleeman BD. Lattice and non-lattice models of tumour angiogenesis. Bull. Math. Biol. 2004;66:1785–819. [PubMed]
419. Plank MJ, Sleeman BD, Jones PF. A mathematical model of tumour angiogenesis, regulated by vascular endothelial growth factor and the angiopoietins. J. Theor. Biol. 2004;229:435–54. [PubMed]
420. Please CP, Pettet GJ, McElwain DLS. A new approach to modeling the formation of necrotic regions in tumors. Appl. Math. Lett. 1998;11:89–94.
421. Please CP, Pettet GJ, McElwain DLS. Avascular tumour dynamics and necrosis. Math. Models Appl. Sci. 1999;9:569–79.
422. Poplawski NJ, Agero U, Gens JS, Swat M, Glazier JA, Anderson ARA. Front instabilities and invasiveness of simulated avascular tumors. Bull. Math. Biol. 2009;71:1189–227. [PMC free article] [PubMed]
423. Postovit LM, Adams MA, Lash GE, Heaton JP, Graham CH. Oxygenmediated regulation of tumor cell invasiveness involvement of a nitric oxide signaling pathway. J. Biol. Chem. 2002;277:35730–7. [PubMed]
424. Pouysségur J, Dayan F, Mazure NM. Hypoxia signalling in cancer and approaches to enforce tumour regression. Nature. 2006;441:437–43. [PubMed]
425. Prall F. Tumour budding in colorectal carcinoma. Histopathology. 2007;50:151–62. [PubMed]
426. Preusser M, Heinzl H, Gelpi E, Schonegger K, Haberler C, Birner P, Marosi C, Hegi M, Gorlia T, Hainfellner JA. Histopathologic assessment of hot-spot microvessel density and vascular patterns in glioblastoma: poor observer agreement limits clinical utility as prognostic factors: a translational research project of the european organization for research and treatment of cancer brain tumor group. Cancer. 2006;107:162–70. [PubMed]
427. Preziosi L. Cancer Modelling and Simulation. Chapman and Hall/CRC; London: 2003.
428. Preziosi L, Ambrosi D, Verdier C. An elasto-visco-plastic model of cell aggregates. J. Theor. Biol. 2009 (doi:10.16/j.jtbi.2009.08.023) [PubMed]
429. Preziosi L, Astanin S. Modelling the formation of capillaries. In: Quarteroni A, editor. Complex Systems in Biomedicine. Springer; Milan: 2006.
430. Preziosi L, Tosin A. Multiphase and multiscale trends in cancer modeling. Math. Model. Nat. Phenom. 2009;4:1–11.
431. Preziosi L, Tosin A. Multiphase modeling of tumor growth and extracellular matrix interaction: mathematical tools and applications. J. Math. Biol. 2009;58:625–56. [PubMed]
432. Pries AR, Reglin B, Secomb TW. Structural adaptation and stability of microvascular networks: functional roles of adaptive responses. Am. J. Physiol. Heart Circ. Physiol. 2001;281:H1015–25. [PubMed]
433. Pries AR, Reglin B, Secomb TW. Structural adaptation of vascular networks: the role of pressure response. Hypertension. 2001;38:1476–9. [PubMed]
434. Pries AR, Secomb TW. Control of blood vessel structure: insights from theoretical models. Am. J. Physiol. Heart Circ. Physiol. 2005;288:1010–15. [PubMed]
435. Pries AR, Secomb TW. Modeling structural adaptation of microcirculation. Microcirculation. 2008;15:753–64. [PMC free article] [PubMed]
436. Pries AR, Secomb TW, Gaehtgens P. Biophysical aspects of blood flow in the microvasculature. Cardiovasc. Res. 1996;32:654–67. [PubMed]
437. Pries AR, Secomb TW, Gaehtgens P. Structural adaptation and stability of microvascular networks: theory and simulations. Am. J. Physiol. Heart Cir. Physiol. 1998;275:H349–60. [PubMed]
438. Quaranta V, Rejniak K, Gerlee P, Anderson ARA. Invasion emerges from cancer cell adaptation to competitive microenvironments: quantitative predictions from multiscale mathematical models. Sem. Cancer Biol. 2008;18:338–48. [PMC free article] [PubMed]
439. Quaranta V, Weaver AM, Cummings PT, Anderson ARA. Mathematical modeling of cancer: the future of prognosis and treatment. Clin. Chim. Acta. 2005;357:173–9. [PubMed]
440. Radszuweit M, Block M, Hengstler JG, Schöll E, Drasdo D. Comparing the growth kinetics of cell populations in two and three dimensions. Phys. Rev. E. 2009;79:051907. [PubMed]
441. Ramanathan A, Wang C, Schreiber SL. Perturbational profiling of a cell-line model of tumorigenesis by using metabolic measurements. Proc. Natl. Acad. Sci. 2005;102:5992–7. [PubMed]
442. Ramis-Conde I, Chaplain MAJ, Anderson ARA. Mathematical modelling of cancer cell invasion of tissue. Math. Comput. Model. 2008;47:533–45.
443. Ramis-Conde I, Chaplain MAJ, Anderson ARA, Drasdo D. Multi-scale modeling of cancer cell intravasation: the role of cadherins in metastasis. Phys. Biol. 2009;6:016008. [PubMed]
444. Ramis-Conde I, Drasdo D, Anderson ARA, Chaplain MAJ. Modeling the inuence of the e-cadherin-beta-catenin pathway in cancer cell invasion: a multiscale approach. Biophys. J. 2008;95:155–65. [PubMed]
445. Rejniak K. A single-cell approach in modeling the dynamics of tumor microregions. Math. Biosci. Eng. 2005;2:643–55. [PubMed]
446. Rejniak K. An immersed boundary framework for modeling the growth of individual cells: an application to the early tumour development. J. Theor. Biol. 2007;247:186–204. [PubMed]
447. Rejniak K, Anderson ARA. A computational study of the development of epithelial acini: I. Sufficient conditions for the formation of a hollow structure. Bull. Math. Biol. 2008;70:677–712. [PMC free article] [PubMed]
448. Rejniak K, Anderson ARA. A computational study of the development of epithelial acini: II. Necessary conditions for structure and lumen stability. Bull. Math. Biol. 2008;70:1450–79. [PMC free article] [PubMed]
449. Rejniak K, Dillon R. A single cell-based model of the ductal tumor microarchitecture. Comput. Math. Methods Med. 2007;8:51–69.
450. Ribba B, Alarcón T, Marron K, Maini PK, Agur Z. The use of hybrid cellular automaton models for improving cancer therapy. In: Chopard B, et al., editors. ACRI, LNCS. Springer; Berlin: 2004. pp. 444–53.
451. Ribba B, Colin T, Schnell S. A multiscale mathematical model of cancer, and its use in analyzing irradiation therapies. Theor. Biol. Med. Modelling. 2006;3:7. [PMC free article] [PubMed]
452. Ridley AJ, Schwartz MA, Burridge K, Firtel RA, Ginsberg MH, Borisy G, Parsons JT, Horwitz AR. Cell migration: Integrating signals from front to back. Science. 2003;302:1704–09. [PubMed]
453. Rini BI. Vegf-targeted therapy in metastatic renal cell carcinoma. Oncologist. 2005;10:191–7. [PubMed]
454. Rofstad E, Halsør E. Hypoxia-associated spontaneous pulmonary metastasis in human melanoma xenographs: involvement of microvascular hotspots induced in hypoxic foci by interleukin. Br. J. Cancer. 2002;86:301–08. [PMC free article] [PubMed]
455. Rofstad E, Rasmussen H, Galappathi K, Mathiesen B, Nilsen K, Graff BA. Hypoxia promotes lymph node metastasis in human melanoma xenografts by up-regulating the urokinase-type plasminogen activator receptor. Cancer Res. 2002;62:1847–53. [PubMed]
456. Roose T, Chapman SJ, Maini PK. Mathematical models of avascular tumor growth. SIAM Rev. 2007;49:179–208.
457. Roose T, Netti PA, Munn LL, Boucher Y, Jain R. Solid stress generated by spheroid growth using a linear poroelastic model. Microvascular Res. 2003;66:204–12. [PubMed]
458. Rubenstein BM, Kaufman LJ. The role of extracellular matrix in glioma invasion: A cellular potts model approach. Biophys. J. 2008;95:5661–80. [PubMed]
459. Rubenstein JL, Kim J, Ozawa T, Zhang K, Westphal M, Deen DF, Shuman MA. Anti-vegf antibody treatment of glioblastoma prolongs survival but results in increased vascular cooption. Neoplasia. 2000;2:306–14. [PMC free article] [PubMed]
460. Rugo HS. Bevacizumab in the treatment of breast cancer: rationale and current data. Oncologist. 2004;1(9 Suppl.):43–9. [PubMed]
461. Sahai E. Mechanisms of cancer cell invasion. Curr. Opin. Genet. Dev. 2005;15:87–96. [PubMed]
462. Sakamoto G. Infiltrating carcinoma: major histological types. In: Page DL, Anderson TJ, editors. Diagnostic Histopathology of the Breast. Churchill-Livingstone; London: 1987.
463. Sakorafas GH, Tsiotou AGH. Ductal carcinoma in situ (dcis) of the breast: evolving perspectives. Cancer Treatment Rev. 2000;26:103–25. [PubMed]
464. Sanga S, Frieboes HB, Zheng X, Gatenby R, Bearer EL, Cristini V. Predictive oncology: a review of multidisciplinary, multiscale in silico modeling linking phenotype, morphology and growth. NeuroImage. 2007;37:S120–34. [PMC free article] [PubMed]
465. Sanga S, Sinek JP, Frieboes HB, Ferrari M, Fruehauf JP, Cristini V. Mathematical modeling of cancer progression and response to chemotherapy. Expert Rev. Anticancer Ther. 2006;6:1361–76. [PubMed]
466. Sansone BC, Del Santo PP, Magnano M, Scalerandi M. Effects of anatomical constraints on tumor growth. Phys. Rev. E. 2002;64:21903.
467. Sansone BC, Scalerandi M, Condat CA. Emergence of taxis and synergy in angiogenesis. Phys. Rev. Lett. 2001;87:128102. [PubMed]
468. Sarntinoranont M, Rooney F, Ferrari M. Interstitial stress and fluid pressure within a growing tumor. Ann. Biomed. Eng. 2003;31:327–35. [PubMed]
469. Savill NJ, Hogeweg P. Modeling morphogenesis: From single cells to crawling slugs. J. Theor. Biol. 1997;184:229–35.
470. Schaller G, Meyer-Hermann M. Kinetic and dynamic delaunay tetrahedralizations in three dimensions. Comput. Phys. Commun. 2004;162:9–23.
471. Schaller G, Meyer-Hermann M. Multicellular tumor spheroid in an off-lattice voronoi-delaunay cell model. Phys. Rev. E. 2005;71:051910. [PubMed]
472. Schmitz JE, Kansal AR, Torquato S. A cellular automaton model of brain tumor treatment and resistance. J. Theor. Med. 2002;4:223–39.
473. Scianna M, Merks R, Preziosi L, Medico E. Individual cell-based models of cell scatter of aro and mlp-29 cells in response to hepatocyte growth factor. J. Theor. Biol. 2009;260:151–60. [PubMed]
474. Seftor EA, Meltzer PS, Kirshmann DA, Pe'er J, Maniotis AJ, Trent JM, Folberg R, Hendrix MJ. Molecular determinants of human uveal melanoma invasion and metastasis. Clin. Exp. Metastasis. 2002;19:233–46. [PubMed]
475. Serini G, Ambrosi D, Giraudo E, Gamba A, Preziosi L, Bussolino F. Modeling the early stages of vascular network assembly. EMBO J. 2003;22:1771–9. [PubMed]
476. Setayeshgar S, Gear CW, Othmer HG, Kevrekidis IG. Application of coarse integration to bacterial chemotaxis SIAM. Multiscale Model. Simul. 2005;4:307–27.
477. Shannon MA, Rubinsky B. The effect of tumor growth on the stress distribution in tissue. Adv. Biol. Heat Mass Transfer. 1992;231:35–8.
478. Sherratt JA. Traveling wave solutions of a mathematical model for tumor encapsulation. SIAM J. Appl. Math. 1999;60:392–407.
479. Sherratt JA. Wavefront propagation in a competition equation with a new motility term modelling contact inhibition between cell populations. Proc. R Soc. 2000;456:2365–86.
480. Sherratt JA, Chaplain MAJ. A new mathematical model for avascular tumour growth. J. Math. Biol. 2001;43:291–312. [PubMed]
481. Shweiki D, Itin A, Soffer D, Keshet E. Vascular endothelial growth factor induced by hypoxia may mediate hypoxia-initiated angiogenesis. Nature. 1992;359:843–45. [PubMed]
482. Shymko RM, Glass L. Cellular and geometric control of tissue growth and mitotic instability. J. Theor. Biol. 1976;63:355–74. [PubMed]
483. Sierra A. Metastases and their microenvironments: linking pathogenesis and therapy. Drug Resist. Updates. 2005;8:247–57. [PubMed]
484. Sinek JP, Frieboes HB, Zheng X, Cristini V. Two-dimensional chemotherapy simulations demonstrate fundamental transport and tumor response limitations involving nanoparticles. Biomed. Microdevices. 2004;6:297–309. [PubMed]
485. Sinek JP, Sanga S, Zheng X, Frieboes HB, Ferrari M, Cristini V. Predicting drug pharmacokinetics and effect in vascularized tumors using computer simulation. J. Math. Biol. 2009;58:485–510. [PMC free article] [PubMed]
486. Smallbone K, Gatenby RA, Gillies RJ, Maini PK, Gavaghan DJ. Metabolic changes during carcinogenesis: Potential impact on invasiveness. J. Theor. Biol. 2007;244:703–13. [PubMed]
487. Smallbone K, Gatenby RA, Maini PK. Mathematical modelling of tumour acidity. J. Theor. Biol. 2008;255:106–12. [PubMed]
488. Smallbone K, Gavaghan DJ, Gatenby RA, Maini PK. The role of acidity in solid tumor growth and invasion. J. Theor. Biol. 2005;235:476–84. [PubMed]
489. Smallbone K, Gavaghan DJ, Maini PK, Brady JM. Quiescence as a mechanism for cyclical hypoxia and acidosis. J. Math. Biol. 2007;55:767–79. [PubMed]
490. Spencer VA, Xu R, Bissell MJ. Extracellular matrix, nuclear and chromatin structure, and gene expression in normal tissues and malignant tumors: a work in progress. Adv. Cancer Res. 2007;97:275–94. [PMC free article] [PubMed]
491. Stamper IJ, Byrne HM, Owen MR, Maini PK. Modelling the role of angiogenesis and vasculogenesis in solid tumour growth. Bull. Math. Biol. 2007;69:2737–72. [PubMed]
492. Steeg PS. Angiogenesis inhibitors: motivators of metastasis? Nature Med. 2003;9:822–3. [PubMed]
493. Stein AM, Demuth T, Mobley D, Berens M, Sander LM. A mathematical model of glioblastoma tumor spheroid invasion in a three-dimensional in vitro experiment. Biophys. J. 2007;92:356–65. [PubMed]
494. Stephanou A, McDougall SR, Anderson ARA, Chaplain MAJ. Mathematical modelling of flow in 2d and 3d vascular networks: applications to anti-angiogenic and chemotherapeutic drug strategies. Math. Comput. Model. 2005;41:1137–56.
495. Stephanou A, McDougall SR, Anderson ARA, Chaplain MAJ. Mathematical modeling of the influence of blood rheological properties upon adaptive tumour-induced angiogenesis. Math. Comput. Model. 2006;44:96–123.
496. Stewart JM, Broadbridge P, Goard JM. Symmetry analysis and numerical modelling of invasion by malignant tumour tissue. Nonlinear Dyn. 2002;28:175–93.
497. Stokes CL, Lauffenburger DA. Analysis of the roles of microvessel endothelial cell random motility and chemotaxis in angiogenesis. J. Theor. Biol. 1991;152:377–403. [PubMed]
498. Stolarska MA, Kim Y, Othmer HG. Multiscale models of cell and tissue dynamics. Phil. Trans. R. Soc. A. 2009;367:3525–53. [PMC free article] [PubMed]
499. Stott EL, Britton N, Glazier JA, Zajac M. Simulation of benign avascular tumour growth using the potts model. Math. Comput. Model. 1999;30:183–98.
500. Stupack DG, Cheresh DA. Get a ligand, get a life: Integrins, signaling and cell survival. J. Cell. Sci. 2002;115:3729–38. [PubMed]
501. Sun CH, Munn LL. Lattice-boltzmann simulation of blood flow in digitized vessel networks. Comput. Math. Appl. 2008;55:1594–600. [PMC free article] [PubMed]
502. Sun S, Wheeler MF, Obeyesekere M, Patrick C., Jr Multiscale angiogenesis modeling using mixed finite element methods. Multiscale Model. Simul. 2005;4:1137–67.
503. Sun S, Wheeler MF, Obeyesekere M, Patrick CW., Jr A deterministic model of growth factor-induced angiogenesis. Bull. Math. Biol. 2005;67:313–37. [PubMed]
504. Sundfor K, Lyng H, Rofstad EK. Tumour hypoxia and vascular density as predictors of metastasis in squamous cell carcinoma of the uterine cervix. Br. J. Cancer. 1998;78:822–7. [PMC free article] [PubMed]
505. Sutherland RM. Cell and environment interactions in tumor microregions: the multicell spheroid model. Science. 1988;240:177–84. [PubMed]
506. Sutherland RM, Carlsson J, Durand RE, Yuhas J. Spheroids in cancer research. Cancer Res. 1981;41:2980–94.
507. Swanson KR, Bridge C, Murray JD, Alvord EC., Jr Virtual and real brain tumors: using mathematical modeling to quantify glioma growth and invasion. J. Neuro. Sci. 2003;216:1–10. [PubMed]
508. Swanson KR, Alvord EC, Jr, Murray JD. A quantitative model for differential motility of gliomas in grey and white matter. Cell Proliferation. 2000;33:317–29. [PubMed]
509. Swanson KR, Alvord EC, Jr, Murray JD. Virtual brain tumours (gliomas) enhance the reality of medical imaging and highlight inadequacies of current therapy. Br. J. Cancer. 2002;86:14–18. [PMC free article] [PubMed]
510. Swanson KR, Alvord EC, Jr, Murray JD. Virtual resection of gliomas: effect of extent of resection on recurrence. J. Math. Comput. Model. 2003;37:1177–90.
511. Swanson KR, Rostomily RC, Alvord EC., Jr A mathematical modeling tool for predicting the survival of individual patients following resection of glioblastoma: a proof of principle. Br. J. Cancer. 2008;98:113–19. [PMC free article] [PubMed]
512. Szeto MD, Chakraborty G, Hadley J, Rockne R, Muzi M, Alvord EC, Jr, Krohn KA, Spence AM, Swanson KR. Quantitative metrics of net proliferation and invasion link biological aggressiveness assessed by mri with hypoxia assessed by fmiso-pet in newly diagnosed glioblastomas. Cancer Res. 2009;69:4502–9. [PMC free article] [PubMed]
513. Szymanska Z, Rodrigo CM, Lachowicz M, Chaplain MAJ. Mathematical modeling of cancer invasion of tissue: the role and effect of nonlocal interactions. Math. Models Methods Appl. Sci. 2009;19:257–81.
514. Szymanska Z, Urbanski J, Marciniak-Czochra A. Mathematical modeling of the influence of heat shock proteins on cancer invasion of tissue. J. Math. Biol. 2009;58:819–44. [PubMed]
515. Tao Y, Chen M. An elliptic-hyperbolic free boundary problem modelling cancer therapy. Nonlinearity. 2006;19:419–40.
516. Tao Y, Guo Q. The competitive dynamics between tumor cells, a replicationcompetent virus and an immune response. J. Math. Biol. 2005;51:37–74. [PubMed]
517. Tao Y, Yoshida N, Guo Q. Nonlinear analysis of a model of vascular tumour growth and treatment. Nonlinearity. 2004;17:867–95.
518. Thomlinson RH, Gray LH. The histological structure of some human lung cancers and the possible implications of radiotherapy. Br. J. Cancer. 1955;9:539–49. [PMC free article] [PubMed]
519. Thorne BC, Bailey AM, Peirce SM. Combining experiments with multi-cell agent-based modeling to study biological tissue patterning. Briefings Bioinformatics. 2007;8:245–57. [PubMed]
520. Tindall MJ, Please CP. Modelling the cell cycle and cell movement in multicellular tumour spheroids. Bull. Math. Biol. 2007;69:1147–65. [PubMed]
521. Tindall MJ, Please CP, Peddie MJ. Modelling the formation of necrotic regions in avascular tumours. Math. Biosci. 2008;211:34–55. [PubMed]
522. Tong S, Yuan F. Numerical simulations of angiogenesis in the cornea. Microvasc. Res. 2001;61:14–27. [PubMed]
523. Tosin A. Multiphase modeling and qualitative analysis of the growth of tumor cords. Networks Heterogen. Media. 2008;3:43–84.
524. Tosin A, Ambrosi D, Preziosi L. Mechanics and chemotaxis in the morphogenesis of vascular networks. Bull. Math. Biol. 2006;68:1819–36. [PubMed]
525. Tracqui P. Biophysical models of tumor growth. Rep. Prog. Phys. 2009;72:056701.
526. Truesdell C, Toupin R. Classical field theories. In: Flugge S, editor. Handbuch der Physik. III/I. Springer; Berlin: 1960.
527. Turner S, Sherratt JA. Intercellular adhesion and cancer invasion: a discrete simulation using the extended potts model. J. Theor. Biol. 2002;216:85–100. [PubMed]
528. Vajkoczy P, Farhadi M, Gaumann A, Heidenreich R, Erber R, Wunder A, Tonn JC, Menger MD, Breier G. Microtumor growth initiates angiogenic sprouting with simultaneous expression of vegf, vegf receptor-2, and angiopoietin-2. J. Clin. Invest. 2002;109:777–85. [PMC free article] [PubMed]
529. van Kempen LCLT, Ruiter DJ, van Muijen GNP, Coussens LM. The tumor microenvironment: a critical determinant of neoplastic evolution. Eur. J. Cell. Biol. 2003;82:539–48. [PubMed]
530. van Leeuwen IMM, Byrne HM, Jensen OE, King JR. Crypt dynamics and colorectal cancer: advances in mathematical modeling. Cell Proliferation. 2006;39:157–81. [PubMed]
531. van Leeuwen IMM, Edwards CM, Ilyas M, Byrne HM. Towards a multiscale model of colorectal cancer. World Gastroenterol. 2007;13:1399–407. [PMC free article] [PubMed]
532. Vukovic V, Haugland HK, Nicklee T, Morrison AJ, Hedley DW. Hypoxia-inducible factor-1 alpha is an intrinsic marker for hypoxia in cervical cancer xenografts. Cancer Res. 2001;61:7394–8. [PubMed]
533. Venkatasubramanian R, Henson MA, Forbes NS. Incorporating energy metabolism into a growth model of multicell spheroids. J. Theor. Biol. 2006;242:440–53. [PubMed]
534. Ventura AC, Jackson TL, Merajver SD. On the role of cell signaling models in cancer research. Cancer Res. 2009;69:400–02. [PMC free article] [PubMed]
535. Visvader JE, Lindeman GJ. Cancer stem cells in solid tumours: accumulating evidence and unresolved questions. Nature Rev. Cancer. 2008;8:755–68. [PubMed]
536. Vollmayr-Lee BP, Rutenberg AD. Stresses in growing soft tissues. Acta Biomater. 2006;2:493–504. [PubMed]
537. Walker C. Global well-posedness of a haptotaxis model with spatial and age structure. Diff. Integ. Eqns. 2007;20:1053–74.
538. Walker C. Global existence for an age and spatially structured haptotaxis model with nonlinear age-boundary conditions. Eur. J. Appl. Math. 2008;19:113–47.
539. Walker C, Webb GF. Global existence of classical solutions for a haptotaxis model. SIAM J. Math. Anal. 2007;38:1694–713.
540. Walles T, Weimer M, Linke K, Michaelis J, Mertsching H. The potential of bioartificial tissues in oncology research and treatment. Onkologie. 2007;30:388–94. [PubMed]
541. Wang Z, Deisboeck TS. Computational modeling of brain tumors: discrete, continuum or hybrid? Scientific Model. Simul. 2008;15:381–93.
542. Wang Z, Zhang L, Sagotsky J, Deisboeck TS. Simulating non-small cell lung cancer with a multiscale agent-based model. Theor. Biol. Med. Model. 2007;4:50. [PMC free article] [PubMed]
543. Ward JP, King JR. Mathematical modelling of avascular tumour growth. IMA J. Math. Appl. Med. Biol. 1997;14:36–69. [PubMed]
544. Ward JP, King JR. Mathematical modelling of avascular-tumour growth: II. modelling growth saturation. Math. Med. Biol. 1999;16:171–211. [PubMed]
545. Ward JP, King JR. Modelling the effect of cell shedding on avasacular tumour growth. J. Theor. Med. 2000;2:155–74.
546. Ward JP, King JR. Mathematical modelling of drug transport in tumour multicell spheroids and monolayer cultures. Math. Biosci. 2003;181:177–207. [PubMed]
547. Wcislo R, Dzwinel W. Particle based model of tumor progression stimulated by the process of angiogenesis. In: Adam J, Bellomo N, editors. Comput. Sci. - ICCS 2008. Springer; Heidelberg: 2008. pp. 177–86.
548. Welter M, Bartha K, Rieger H. Emergent vascular network inhomogenities and resulting blood flow patterns in a growing tumor. J. Theor. Biol. 2008;250:257–80. [PubMed]
549. Welter M, Bartha K, Rieger H. Hot spot formation in tumor vasculature during tumor growth in an arterio-venous-network environment. 2008. arXiv:0801.0654v2[q-bio.To]
550. Wertheimer E, Sasson S, Cerasi E, Ben-Neriah Y. The ubiquitous glucose transporter glut-1 belongs to the glucose-regulated protein family of stress-inducible proteins. Proc. Natl Acad. Sci. 2008;88:2525–9. [PubMed]
551. Whitaker N, Harrington HA, Maier M, Naidoo L, Kevrekidis PG. A hybrid model for tumor-induced angiogenesis in the cornea in the presence of inhibitors. Math. Comput. Modelling. 2007;46:513–24.
552. Wicha MS, Liu S, Dontu G. Cancer stem cells: An old idea- a paradigm shift. Cancer Res. 2006;66:1883–90. [PubMed]
553. Wise SM, Kim JS, Lowengrub JS. Solving the regularized, strongly anisotropic Chan-Hilliard equation by an adaptive nonlinear multigrid method. J. Comput. Phys. 2007;226:414–46.
554. Wise SM, Lowengrub JS, Cristini V. An adaptive algorithm for simulating solid tumor growth using mixture models. J. Math. Comp. Model. submitted. [PMC free article] [PubMed]
555. Wise SM, Lowengrub JS, Frieboes HB, Cristini V. Three-dimensional multispecies nonlinear tumor growth: I. Model and numerical method. J. Theor. Biol. 2008;253:524–43. [PMC free article] [PubMed]
556. Wolf K, Friedl P. Molecular mechanisms of cancer cell invasion and plasticity. Br. J. Dermatol. 2006;154:11–15. [PubMed]
557. Woodward DE, Cook J, Tracqui P, Cruywagen GC, Murray JD, Alvord EC., Jr A mathematical model of glioma growth: The effect of surgical resection. Cell Proliferation. 1996;29:269–88. [PubMed]
558. Wu J, Cui S. Asymptotic behavior of solutions of a free boundary problem modeling the growth of tumors with stokes equations Discrete. Contin. Dyn. Syst. 2009;24:625–51.
559. Wu J, Zhou F, Cui S. Simulation of microcirculation in solid tumors. Complex Med. Eng. CME 2007. IEEE/ICME Int. Conf. (Beijing, China) 2007:1555–63.
560. Wurzel M, Schaller C, Simon M, Deutsch A. Cancer cell invasion of brain tissue: guided by a prepattern? J. Theor. Med. 2005;6:21–31.
561. Xu R, Spencer VA, Bissell MJ. Extracellular matrix-regulated gene expression requires cooperation of swi/snf and transcription factors. J. Biol. Chem. 2007;282:14992–9. [PMC free article] [PubMed]
562. Xu S. Hopf bifurcation of a free boundary problem modeling tumor growth with two time delays. Chaos Solitons Fractals. 2009;41:2491–4.
563. Xu Y. A free boundary problem model of ductal carcinoma in situ 2004. Discrete Contin. Dyn. Syst. 2004;4:337–48.
564. Xu Y, Gilbert R. Some inverse problems raised from a mathematical model of ductal carcinoma in situ 2009. Math. Comput. Model. 2009;49:814–28.
565. Yamaguchi H, Wyckoff J, Condeelis J. Cell migration in tumors. Curr. Opin. Cell Biol. 2005;17:559–64. [PubMed]
566. Yin SP, et al. Maspin retards cell detachment via a novel interaction with the urokinase-type plasminogen activator/urokinase-type plasminogen activator receptor system. Cancer Res. 2006;66:4173–81. [PubMed]
567. Young SD, Hill RP. Effects of reoxygenation of cells from hypoxic regions of solid tumors: anticancer drug sensitivity and metastatic potential. J. Natl. Cancer Inst. 1990;82:338–9. [PubMed]
568. Young SD, Marshall RS, Hill RP. Hypoxia induces dna overreplication and enhances metastatic potential of murine tumor cells. Proc. Natl Acad. Sci. USA. 1988;85:9533–7. [PubMed]
569. Zagzag D, Amirnovin R, Greco MA, Yee H, Holash J, Wiegand SJ, Zabski S, Yancopoulos GD, Grumet M. Vascular apoptosis and involution in gliomas precede neovascularization: a novel concept for glioma growth and angiogenesis. Lab. Invest. 2000;80:837–49. [PubMed]
570. Zajac M, Jones GL, Glazier JA. Model of convergent extension in animal morphogenesis. Phys. Rev. Lett. 2000;85:2022–25. [PubMed]
571. Zhang L, Athale CA, Deisboeck TS. Development of a three-dimensional multiscale agent-based tumor model: simulating gene-protein interaction profiles, cell phenotypes and multicellular patterns in brain cancer. J. Theor. Biol. 2007;244:96–107. [PubMed]
572. Zhang L, Strouthos CG, Wang Z, Deisboeck TS. Simulating brain tumor heterogeneity with a multiscale agent-based model: Linking molecular signatures, phenotypes and expansion rate. Math. Comput. Model. 2009;49:307–19. [PMC free article] [PubMed]
573. Zhang L, Wang Z, Sagotsky JA, Deisboeck TS. Multiscale agent-based cancer modeling. J. Math. Biol. 2009;58:545–59. [PubMed]
574. Zhao G, Wu J, Xu S, Collins MW, Long Q, Koenig CS, Jiang Y, Wang J, Padhani AR. Numerical simulation of blood flow and interstitial fluid pressure in solid tumor microcirculation based on tumor-induced angiogenesis. Mech. Sin. 2007;23:477–83.
575. Zheng X, Wise SM, Cristini V. Nonlinear simulation of tumor necrosis, neovascularization and tissue invasion via an adaptive finite-element/level-set method. Bull. Math. Biol. 2005;67:211–59. [PubMed]
576. Zhou F, Cui S. Bifurcation for a free boundary problem modeling the growth of multi-layer tumors. Nonlinear Anal. Theory Methods Appl. 2008;68:2128–45.
577. Zhou F, Escher J, Cui S. Well-posedness and stability of a free boundary problem modeling the growth of multi-layer tumors. J. Diff. Eqns. 2008;244:2909–33.
578. Zhu H, Yuan W, Ou C. Justification for wavefront propagation in a tumour growth model with contact inhibition. Proc. R. Soc. 2008;464:1257–73.