## Abstract

The hallmark of malignant tumours is their spread into neighbouring tissue and metastasis to distant organs, which can lead to life threatening consequences. One of the defining characteristics of aggressive tumours is an unstable morphology, including the formation of invasive fingers and protrusions observed both *in vitro* and *in vivo*. In spite of extensive biological, clinical and modelling study and research at physical scales ranging from the molecular to the tissue, the driving dynamics of tumour invasiveness are not completely understood, partly because it is challenging to observe and study cancer as a multi-scale system. Mathematical modelling has been applied to provide further insights into these complex invasive and metastatic behaviours. Modelling a solid tumour as an incompressible fluid, we consider three possible constitutive relations to describe tumour growth, namely Darcy's law, Stokes' law and the combined Darcy–Stokes law. We study the tumour morphological stability described by each model and evaluate the consistency between theoretical model predictions and experimental data from *in vitro* three-dimensional multicellular tumour spheroids. The analysis reveals that the Stokes model is the most consistent with the experimental observations, and that it predicts our experimental tumour growth is marginally stable. We further show that it is feasible to extract parameter values from a limited set of data and create a self-consistent modelling framework that can be extended to the multi-scale study of cancer.

## 1. Introduction

We study tumour morphological stability by employing three mathematical models to gain insight into tumour invasion and metastasis, and to evaluate the consistency between the model predictions and experimental observations. There has been extensive progress in employing mathematical modelling to study solid tumour growth over the past few decades, and this work has provided insight into the understanding of experimental and clinical data. Most models fall into two categories, discrete cell-based and continuum models (see recent reviews Byrne *et al.* 2006; Fasano *et al.* 2006; Galle *et al.* 2006; Drasdo & Hohme 2007; Friedman *et al.* 2007; Graziano & Preziosi 2007; van Leeuwen *et al.* 2007; Martinsa *et al.* 2007; Roose *et al.* 2007; Sanga *et al.* 2007; Anderson & Quaranta 2008; Bellomo *et al.* 2008; Wang & Deisboeck 2008; Deisboeck *et al.* 2009; Preziosi & Tosin 2009; Tracqui 2009; Ventura *et al.* 2009; Zhang *et al.* 2009; Lowengrub *et al.* 2010). For example, models have been applied to brain cancer (e.g. Swanson *et al.* 2003; Jbabdi *et al.* 2005; Khain & Sander 2006) and breast cancer (e.g. Franks *et al.* 2003*a*,*b*). Furthermore, modelling has shown that tumour morphology may serve as a predictor of invasiveness (Cristini *et al.* 2003, 2005; Frieboes *et al.* 2006; Friedman & Hu 2007; Macklin & Lowengrub 2007; Anderson & Quaranta 2008; Lowengrub *et al.* 2010). In particular, Byrne & Chaplain (1996, 1997) revealed that cell–cell adhesion and external nutrient concentration are key parameters controlling the stability of three-dimensional multicellular spheroids. Extending this work, Cristini *et al.* (2003) proposed that two non-dimensional parameters could help describe the complexity of tumour shapes observed in both avascular (e.g. *in vitro*) and vascular (e.g. *in vivo*) conditions. Here, we complement and extend the work in Cristini *et al.* (2003) by employing three macroscale single-phase sharp interface models, and by applying linear stability analyses to further elucidate the determinants of tumour shape.

Greenspan (1976) considered necrotic tumours in the avascular stage, where growth is regulated solely by nutrient in the surrounding micro-environment. Byrne & Chaplain (1995) proposed a model for non-necrotic tumours where nutrient is supplied through the surrounding vascularized environment. Following this work (Byrne & Chaplain 1995) and our previous work (Cristini *et al.* 2003), we first consider non-necrotic tumours in a vascularized environment, and later simplify the model to the avascular condition to compare its predictions to our *in vitro* experiments (Frieboes *et al.* 2006). We employ Darcy's law and Stokes flow as constitutive laws in describing the deformation and stress fields of the tissue. Combining Darcy's law with Stokes flow gives a third constitutive relation, Darcy–Stokes flow (also known as the Brinkman equation, Zheng *et al.* 2005). We describe cell–cell adhesive forces by a surface tension at the tumour–tissue interface. Tumour growth is governed by the balance between cell mitosis and apoptosis, as well as cell–cell adhesion. The rate of mitosis depends on the concentration of nutrient that obeys a diffusion–reaction equation within the tumour volume.

The Darcy model, which models flow through a porous medium, was previously considered by Greenspan (1976), Byrne & Chaplain (1995), Friedman & Reitich (1999), Cristini *et al.* (2003), and others, while Stokes' law was studied by Friedman & Hu (2007). Both models were investigated by King & Franks (2004, 2007*a*,*b*) and Franks & King (2003). The combined Darcy–Stokes law was simulated but not analysed in Zheng *et al.* (2005). Cristini *et al.* (2003) showed that a non-dimensionalization of the Darcy model gives rise to two dimensionless parameters and , respectively, representing the ratio between apoptosis and proliferation, and the relative rate of proliferation to the relaxation mechanisms (cell mobility and cell–cell adhesion). Linear stability analysis reveals that these are competing forces that control tumour morphology. Frieboes *et al.* (2006) applied these results to the study of glioblastoma tumour spheroids *in vitro* using glucose and growth serum as moderating factors of cell adhesion and cell proliferation. They found an estimated experimental range of and within which the tumour morphology was predicted marginally stable, in excellent agreement with the experimental observations. In addition, Cristini *et al.* (2003) showed that the parameter and another dimensionless parameter *B*, representing the effect of vascularization, subdivide tumour growth into regimes of low, moderate and high vascularization. In particular, there exist critical conditions in the low vascularization stage under which unstable tumour morphologies may develop. Under isotropic vascularization, growth in the moderate and high vascularization regimes is always stable, indicating that additional inhomogeneities are needed for unstable growth. Here, we perform an analogous analysis using the Stokes and Darcy–Stokes models.

Unlike the Darcy model, the Stokes and Darcy–Stokes models account for viscous stress. We assess analytical results against the experiments in Frieboes *et al.* (2006) to determine the consistency between the model predictions of shape stability and the experimental data. Following Cristini *et al.* (2003) and using the theory of vector spherical harmonics (Hill 1954), we perform a linear stability analysis for each of the three models. The results give rise to a marginally stable value _{M}, which depends on the tumour radius and wavenumber (and in the Darcy and Darcy–Stokes models). This value divides the tumour growth into two regions describing stable and unstable morphologies, thus defining the marginal stability curve. Comparison between model predictions of the value of for the experiments, which we term _{P}, and _{M} reveals whether the model predicts stable or unstable growth. Deviations of _{p} from _{M} are quantified to give analytical stability measures, which can be assessed against experimental measurements to determine which model is more consistent, hence more predictive of tumour morphology.

The outline of this paper is as follows. In §2, we introduce the mathematical models and present the corresponding results employed to assess the models. We also describe the process of the qualitative assessment of the *in vitro* tumour spheroid cultures. In §3, we present the analysis of shape stability predicted by each model and the assessment of each one against the experimental results. We end with conclusions and discussion of future work in §4.

## 2. Methods

We briefly discuss the experiments used to obtain observations of *in vitro* tumour morphologies. We then introduce and analyse three single-phase sharp interface models to compare their morphology stability predictions with experimental observations of *in vitro* tumour spheroids. In the models, we assume that the tumour cell population is homogeneous. We treat the tumour as an incompressible fluid and employ different constitutive laws to model tissue stresses. Growth is regulated by cell substrates, e.g. oxygen and glucose, which regulate biophysical processes such as cell proliferation and apoptosis. Cell–cell adhesion helps to maintain tumour compactness, and pressure owing to cell proliferation acts as an expansive force.

### 2.1. Cell culture

Experimental data were obtained as part of a study of *in vitro* and *in silico* tumour invasion, as detailed in Frieboes *et al.* (2006). Briefly, the ACBT (grade IV human glioblastoma multiforme) cell line was cultured to grow tumour spheroids using a liquid-overlay technique (Sutherland *et al.* 1981). Three levels of each glucose (4.5, 2.75 and 1.0 g l^{−1}) and foetal bovine serum (FBS; 10%, 5% and 1%) were used to obtain nine sets of measurements. Each set consisted of 48 tumour growth measurements. Glucose mainly affects cell adhesion and metabolism, while serum mainly affects proliferation. In tumour spheroids greater than 1 mm in diameter, cell adhesion may decrease with glucose concentration because higher levels of glucose lead to reduced oxygenation (Mueller-Klieser *et al.* 1986), thus increasing tumour cell motility as a response to hypoxia (Postovit *et al.* 2002).

Spheroid growth was observed using a Leitz microscope at magnification 100× over a period of 36 days. Spheroid diameter was measured via the arithmetic average of two orthogonal parameters, one being the longest observable axis. Photographs were taken with an Olympus camera with a photography window of 1130 by 1430 µm.

### 2.2. Governing equations

Following Byrne & Chaplain (1995), Friedman & Reitich (1999) and Cristini *et al.* (2003) and others, we present the governing equations for growth for a non-necrotic tumour in a vascularized environment. Let *Ω* = *Ω*(*t*) be the tumour volume and *Σ* = *Σ*(*t*) be the sharp interface separating the tumour from the host (or cell culture medium). We assume that the concentration of cell substrates (e.g. oxygen and nutrients such as glucose and growth factors) *σ*(**x**, *t*) in *Ω* satisfies the reaction–diffusion equation:
2.1
where is the diffusion coefficient and *Γ* is the rate at which cell substrates are added to *Ω* accounting for all sources and sinks of substrates in the tumour volume.

These substrates are supplied by a vasculature at a rate *Γ*_{B} (*σ*, *σ*_{B}), where *σ*_{B} is the concentration (uniform) in the blood:
2.2
Hence *Γ* is given by
2.3
where *λ*_{σ} is the uptake rate. We assume constant and uniform cell substrate concentration at the interface *Σ* (and outside *Ω*),
2.4

As in Greenspan (1976), Byrne & Chaplain (1995), Friedman & Reitich (1999) and Cristini *et al.* (2003) and others, we assume constant cell density in the tumour domain. Hence, mass changes correspond to volume changes. Defining **v** to be the cell velocity, the local rate of volume change ∇ · **v** is given by
2.5
where *λ*_{P} is the net cell-proliferation rate depending on the rate of mitosis *b* and a measure of apoptosis *λ*_{A}, both assumed to be uniform:
2.6

### 2.3. The Darcy model

We consider different constitutive laws in describing the deformation and stress of the tissue to model the tumour as various types of fluid.

As an ideal fluid, we use Darcy's law (Cristini *et al.* 2003):
2.7
where *α*^{−1} is the cell motility. At the boundary, the pressure is assumed to satisfy the Laplace–Young boundary condition (Young 1805),
2.8
where *τ* is the surface tension used to model cell–cell adhesion and *κ* is the local total curvature. This boundary condition is used in fluid mechanics to model force imbalance between two heterogeneous phases. We assume here that tumour and host cells tend to stay with their own kind, causing an imbalance of adhesion forces at the interface.

Finally, we assume the normal velocity at the interface *Σ* to be
2.9

The tumour interface *Σ* evolves according to the normal velocity:
2.10

Applying the theory of vector spherical harmonics (Hill 1954) to the linear analysis of each model, we first recapitulate the linear results for the Darcy model in Cristini *et al.* (2003). We also use a different non-dimensionalization than that in Cristini *et al.* (2003).

#### 2.3.1. Non-dimensionalization

We define the following non-dimensional parameters
2.11
where *L* is the length scale, *λ*_{t} is the time scale, *σ*_{∞} is the cell substrate characteristic, and *P*_{a} is the pressure characteristic value. We define the effective mitosis rate to be proportional to the cell substrate level at the interface, i.e. *λ*_{M} = *b**σ*_{∞}, and from equation (2.5) we obtain the intrinsic time scale, *λ*_{t} = *λ*_{M}. From equations (2.1)–(2.3), the intrinsic diffusional length scale is *L* = ^{1/2} (*λ*_{B} + *λ*_{σ})^{−1/2}. Finally, the pressure scale can be taken from Darcy's law as *P*_{a} = *α**L*^{2}*λ*_{M} (equation (2.7)).

We assume *λ*_{M} *L*^{2}/≪1, i.e. nutrient diffusion is much faster than mitosis. This assumption allows an approximation of the dimensional reaction–diffusion equation in equation (2.1) by a quasi-steady condition in equation (2.16).

We introduce the modified concentration *σ̄* such that
2.12
where *B* is the dimensionless relative effect of vascularization,
2.13

Dropping the bars and primes, as in Cristini *et al.* (2003), we obtain the non-dimensional Darcy system:
2.14
2.15
2.16
with boundary conditions
2.17
2.18
2.19

The interface evolution is as in equation (2.10). This system introduces two new dimensionless parameters and , representing the non-dimensional relative rate of apoptosis and relative strength of cell–cell adhesion, respectively: 2.20

We now summarize the linear stability analysis in Cristini *et al.* (2003) and Li *et al.* (2007). From equation (2.10), which assumes that the tumour interface evolves in the direction of the normal velocity given in equation (2.9), the evolution equation for the unperturbed tumour radius *r*(*t*) = *R* is
2.21
which depends on the apoptosis parameter and vasculature parameter *B*. The right-hand side is the radial component of the velocity field, which can be obtained from the Darcy system.

We perturb the spherical tumour interface *Σ* as follows:
2.22
where *δ* is the dimensionless perturbation size, *Y*_{l,m} a spherical harmonic, *l* and *θ* the polar wavenumber and angle, and *m* and *ϕ* the azimuthal wavenumber and angle. As in Cristini *et al.* (2003) and Li *et al.* (2007), the evolution equation for the shape perturbation *δ*/*R* of the interface is
2.23
where *I*_{l} (*R*) is the modified Bessel function of the first kind. Note that *δ*/*R* rather than *δ* is appropriate to characterize the perturbations because the tumour radius *R* varies in time. Observe that the right-hand side of the equation increases with increasing (high cell death) and decreases with increasing ^{−1} (high cell adhesion), implying that promotes shape instability while promotes stable morphology. This suggests that morphology is determined by the competition between cell death and cell–cell adhesion.

### 2.4. Stokes model

Modelling the tumour as a viscous fluid, we use Stokes flow to describe the cell velocity (Franks & King 2003; King & Franks 2007*a*),
2.24
where **T** is the stress tensor defined as
2.25

Here, **T** takes into account the local rate of strain, dilatation and pressure. The parameters *ν*_{1} and *ν*_{2} are the viscosity coefficients and *p* is the pressure. As in the Darcy model, the stress tensor is assumed to satisfy the Laplace–Young boundary condition at the tumour interface:
2.26

The normal velocity at the interface is given by equation (2.9), and the tumour interface evolves as in equation (2.10).

#### 2.4.1. Non-dimensionalization

We rescale the stress tensor by a characteristic value *T*_{a}:
2.27

We use the same non-dimensional scaling for the other variables as in the Darcy model. From equations (2.24) and (2.25), we take *T*_{a} = *λ*_{M}*ν*_{1}, where *λ*_{M} is the mitosis rate. We define the following non-dimensional parameters,
2.28
where is the relative strength of cell–cell adhesion as before and *ν* is the ratio between the two viscosity coefficients. Note that the definition of here is different from that of the Darcy model, equation (2.20).

Summarizing the non-dimensional Stokes system, we have
2.29
2.30
2.31
with boundary conditions
2.32
2.33
2.34
where **T** = ∇**v** + (∇**v**)^{T} + *ν*(∇ · **v**)**I** − *p***I** is the non-dimensional stress tensor. The interface evolves according to the normal velocity as in equation (2.10).

The non-dimensional Stokes system introduces a new dimensionless parameter *ν* in addition to and , where is the same as in the Darcy model.

#### 2.4.2. Linear stability analysis

We perform a linear stability analysis of the non-dimensional Stokes system. For an unperturbed spherical tumour interface, we obtain the same evolution equation for the tumour radius as before in equation (2.21). For the perturbed spherical tumour interface, we obtain a new evolution equation:
2.35
where *H*(*R*) = (1/tanh *R*) − (1/*R*). The functions *I*_{l}(*R*) and *j*_{l}(i*r*) (where ) are the modified Bessel function of the first kind and the spherical Bessel function, respectively, related by *j*_{l} (i*x*) = *e*^{(1/2)lπi}*I*_{l+1/2}(*x*). Note that the spherical Bessel function is a solution of the radial part of the nutrient equation, namely equation (2.31), which is a form of the Helmholtz equation (Abramowitz & Stegun 1972). The term *j*″_{l} (i*R*)/*j*_{l} (i*R*) is real and can be written in terms of *I*_{l} (*R*) as
2.36

As in the Darcy model, cell adhesion (^{−1}) promotes a stable morphology. However, the shape evolution equation (2.35) is independent of , which means that cell apoptosis does not directly affect shape morphology. Indirectly, the apoptosis parameter () may affect shape morphology through equation (2.21) for the tumour radius, suggesting that it may play a role in promoting shape instability. Our results below show that the marginal stability curve _{M} (*l,,R*) is insensitive to *l* and *λ*. Thus, the main controlling factors of shape morphology in the Stokes model are and the tumour radius *R*.

### 2.5. The Darcy–Stokes model

We now combine the Darcy and Stokes models by adding a dissipation or friction term to Stokes' law (or equivalently a viscous stress to Darcy's law):
2.37
where **T** is the viscous stress tensor defined as in the Stokes model. The constant *α*^{−1} > 0 represents the cell motility as in the Darcy model. Equation (2.37) is also known as the Brinkman equation, considered by Zheng *et al.* (2005). We impose the same Laplace–Young boundary condition as in the Stokes model,
2.38

As before, we assume the same normal velocity at the tumour boundary 2.39 and the tumour interface evolution as in equation (2.10).

#### 2.5.1. Non-dimensionalization

We can non-dimensionalize the Darcy–Stokes model using the dimensionless parameters in either the Darcy or the Stokes model. Depending on the choice of the parameter set, we find that the Darcy–Stokes model either converges to the Stokes model as *α* → 0 or to the Darcy model as *α* → ∞.

Using the parameter set as in the Stokes model, we arrive at the following non-dimensional system: 2.40 2.41 2.42 with boundary conditions 2.43 2.44 2.45

As before, **T** is the non-dimensional viscous stress tensor. The tumour interface evolution is as in equation (2.10). Here, *α̃*^{−1} = *ν*_{1}/*α**L*^{2} is a new dimensionless parameter representing cell motility. This system differs from the non-dimensional Stokes system by the non-dimensional Brinkman equation (2.41).

#### 2.5.2. Linear stability analysis

We perform the linear stability analysis of the Darcy–Stokes model. By solving the non-dimensional system on an unperturbed spherical interface, we obtain the same evolution equation for the unperturbed tumour radius *R* as in the previous two models (equation (2.21)).

In the presence of a perturbed interface, the evolution equation for the shape perturbation *δ*/*R* is
2.46
where
2.47
2.48
2.49

The shape evolution depends on both parameters and . In particular, the dependence on arises from the Darcy component of the model. Here we replaced *α̃* with *α*. Note that the rational functions with the imaginary number *i* involving spherical Bessel functions can be converted into real rational functions involving modified Bessel functions of the first kind.

## 3. Results

### 3.1. Stability analysis of avascular growth

In this section we study the stability of tumour growth in the avascular stage, i.e. *B* = 0.

#### 3.1.1. The Darcy model

A marginally stable (or critical) value of the adhesion parameter _{M,D} (*l,,R*) is obtained by setting the time derivative of the shape factor equation (2.23) to zero. We consider first the case = 0 and plot _{M,D} versus the unperturbed tumour radius *R* for different modes *l* in figure 1*a*. For a single mode *l*, the _{M,D} curve separates the plane into upper and lower regions describing stable and unstable shape morphologies, respectively. The curves lie below the line _{M,D} = 0, implying that shape evolution is always stable since is non-negative. Hence, growth remains stable for = 0. As indicated by equation (2.23), promotes shape instability. Considering instead critical values such that d*R*/d*t* = 0 from equation (2.21), we plot _{M,D} versus the unperturbed tumour radius *R* for different modes *l* in figure 1*b*. The curves now lie above the line _{M,D} = 0. Unlike the case = 0, the region above the curve describes an unstable morphology while the region below the curve describes a stable morphology. The curves shift upward with increasing *l*, implying that the model is more stable as *l* increases. This can be seen from equation (2.23), since for fixed and *R*, (*δ*/*R*)^{−1} (d/d*t*)(*δ*/*R*) ∼ −^{−1} *l*^{3} for large *l*, implying that the shape perturbation is significantly more stable as *l* increases.

#### 3.1.2. The Stokes model

We present an analogous analysis for the Stokes model. A marginally stable (or critical) value of _{M,S}(*λ*,*l,R*) for the Stokes model is obtained by setting (d/d*t*)(*δ*/*R*) = 0 in equation (2.35). In figure 2, we plot _{M,S} versus the unperturbed tumour radius *R* for different *l*. Parameter tests (not shown) show that _{M,S} curves are insensitive to varying *λ*; thus we take *λ* = 1. For a given mode *l*, the _{M,S} curve divides the parameter space into regions of stable and unstable growth. Unlike the Darcy case, _{M,S} here is independent of the relative strength of apoptosis, . Recall in the Darcy case that growth is stable if = 0. In the Stokes case, this demonstrates that an unstable morphology is possible even if = 0. Observe also that curves for larger modes cross those for smaller modes. The fact that curves for larger modes lie above the low modes for small *R* implies that in the early stage of development larger modes are more stable than smaller modes, but the reverse is true as the tumour grows to larger radii. Nevertheless, the curves appear clustered together, indicating that tumour stability is relatively insensitive to *l* (compared with the Darcy model). This can be seen from equation (2.35), where (*δ*/*R*)^{−1} (d/d*t*)(*δ*/*R*) ∼ −^{−1} *l* for large *l* and fixed *R* in contrast to the *l*^{3} behaviour in the Darcy model. The viscous effects in the Stokes model damp out the dependence of the perturbation amplitudes on the wavenumber *l* compared to the Darcy model.

#### 3.1.3. The Darcy–Stokes model

Next, we study the marginally stable (critical) value _{M,DS} for the Darcy–Stokes model. Figure 3 shows that, for a single mode (*l* = 3), the _{M,DS} curves for different *α* lie between those of the other two models. For simplicity, we consider the case of *α* = 1. Again, *λ* does not affect _{M,DS}, and so we pick *λ* = 1. Similar to the Darcy model, _{M,DS} depends explicitly on (equation (2.46)). In figure 4*a*, we plot _{M,DS} versus the unperturbed tumour radius *R* for the case = 0. The curves for *l* = 2, 3, 4 lie below the line = 0, suggesting that perturbations involving these modes are always stable. On the other hand, the _{M} curves for *l* ≥ 5 have discontinuities, i.e. _{M,DS}^{−1} vanishes at certain points as seen in figure 4*b*. These curves initially lie above the line = 0, but they become negative after the points of discontinuity, suggesting that instability is possible before these points. That is, shape instability is possible for = 0, which was observed previously in the Stokes model, but only for sufficiently large modes and sufficiently small radii. In the case where is chosen as the critical (steady-state) value in equation (2.21), figure 4*c* shows that the curves behave like those in the Darcy model.

### 3.2. Stability analysis of vascular growth

Next we consider vascular growth. As in Cristini *et al.* (2003), the dimensionless parameters and *B* subdivide growth into three regimes of vascularization. We evaluate tumour stability predictions by the models in these regimes.

#### 3.2.1. Low vascularization

> 0 (i.e. *λ*_{A}/*λ*_{M} > *B*) and *B* < 1.

In this case, we can rescale time by 1 − *B* > 0 and reduce the models to the case of *B* = 0 (avascular stage) with redefined parameters. Apoptosis in the reduced models is redefined as and cell adhesion as , which remains non-negative. Thus, the marginally stable _{M} curves for all three modified models are positive as in §3.1, where they divide growth into stable and unstable regions. Hence, all three models predict both stable and unstable evolution in this regime.

#### 3.2.2. Moderate vascularization

≤ 0 (i.e. 0 ≤ *λ*_{A}/*λ*_{M} ≤ *B*) and *B* < 1.

As in the low vascularization regime, we rescale time and modify and by 1 − *B* > 0 as above. Again, the models reduce to the case of *B* = 0 (avascular stage) with redefined parameters. The Darcy and Stokes models predict different results here. The evolution is stable using the Darcy model, but the Stokes model predicts both stable and unstable growth. On other hand, the Darcy–Stokes model predictions vary with *l*, i.e. stable evolution if 2 ≤ *l* ≤ 4 and both stable and unstable evolution if *l* > 4. From the evolution equations for the unperturbed tumour radius *R*, we see that d*R*/d*t* → −*AR*/3 as *R* → ∞, implying that growth is exponential for < 0 and linear for = 0.

#### 3.2.3. High vascularization

*B* > 1.

When *B* > 1, the marginally stable _{M} curves of all the models are negative, meaning that evolution is always stable independent of the models. Hence, unstable growth in this regime would require tumour microenvironment inhomogeneity (e.g. Cristini *et al.* 2005; Frieboes *et al.* 2007, 2010; Macklin & Lowengrub 2007; Bearer *et al.* 2009), which is not modelled here.

### 3.3. Assessment of the models

We next assess the models against experimental observations of *in vitro* tumour growth. Thus, we evaluate theoretical predictions of shape stability in the avascular growth regime (*B* = 0), for which morphology is predicted to be stable and unstable by all the models. In Frieboes *et al.* (2006), glioblastoma spheroids were cultured using three levels of glucose and three levels of growth serum. Glucose mainly affects cell adhesion and metabolism, while serum mainly affects cell proliferation. For each set of glucose and growth serum concentrations, spheroids were allowed to grow to a steady state (see sample tumour photographs in figure 5). We made qualitative stability evaluations of the morphologies (figure 6) based on different criteria, including shape compactness and the amount of cell shedding. The figure shows the average measured values for tumour radii. Median values are also shown, indicating that most tumours in each category had radii that were close to the average. Using the average values as inputs, we quantify model predictions of shape stability by estimating the values of for the experimental conditions and comparing the estimated = _{p} with the marginally stable (or critical) _{M} values of the shape evolution equations obtained in §3.1 with *B* = 0. We thus compare and contrast the analytical and experimental stability measures to establish the consistency between the models and the experimental data.

The analysis in Frieboes *et al.* (2006), also based on the Darcy model, mainly evaluated the effect of apoptosis on the *in vitro* tumour morphological stability. However, the evolution equation for the shape perturbation does not always depend on the rate of apoptosis, as observed in the Stokes model (equation 2.35). Hence, we do not rely on the analysis of apoptosis, and instead study the morphologies using the adhesion forces, represented by the parameter .

#### 3.3.1. Estimation of apoptosis

We consider when the tumour is at steady state, i.e. when d*R*/d*t* = 0, corresponding to the *in vitro* glioblastoma cultures when the tumour spheroids have reached a diffusion-limited steady-state size. That is, we obtain from equation (2.21) with *B* = 0:
3.1
where *R* is the steady-state experimental radius scaled by the empirically observed oxygen diffusion length *L* ≈ 125 µm. From the experiments (Frieboes *et al.* 2006), it is estimated that 6.240 ≤ *R* ≤ 7.528 (tables 1⇔–3), which yields 0.3456 ≤ ≤ 0.4037 using equation (3.1).

Considering the exponential growth regime for small tumours, i.e. approximating d*R*/d*t* near *R* ≈ 0, we obtain the net mitosis rate:
3.2

Then the evolution of the unperturbed dimensional tumour radius *R̄* becomes (d*R̄*/d*t̄*) = *λ̃*_{M}*R̄*, where *R̄* = *RL*.

#### 3.3.2. Estimation of cell adhesion

We next provide an estimate of from the models, which we use for evaluation against the experimental results.

### Darcy model.

From Darcy's law and the Laplace–Young boundary condition, we have 3.3 where the left is the characteristic pressure at the proliferating rim from the dimensional Darcy's law (equation (2.7)) and the right is the characteristic pressure at the tumour boundary from the boundary condition (equation (2.8)).

Equating the two pressures and using the definition of (equation (2.20)), we obtain 3.4

Taking this as the prediction for the adhesion parameter = _{p}, we get
3.5
where *R* = *R̄*/*L* is the non-dimensional tumour radius at steady-state, and is the apoptosis parameter at steady state from equation (3.1).

In Frieboes *et al.* (2006), the effect of was assumed small and neglected. However, we consider it here, since based on the experimental data, 0.34 < < 0.41.

### Stokes model.

From the Stokes equation (2.24) and the corresponding Laplace–Young boundary condition (2.26), we have
3.6
where *T*_{a} are the characteristic stresses in the proliferating rim and tumour boundary. Equating the two and using the definition of ^{−1} in equation (2.28), we obtain
3.7

As before, this provides a prediction for the adhesion parameter = _{p}:
3.8

Note that the prediction for the cell adhesion is the same for both models even though is defined differently. The same prediction is also obtained for the Darcy–Stokes model, showing that this formula is independent of the constitutive model.

#### 3.3.3. Comparison with experiments

Following Frieboes *et al.* (2006), we set the mitosis rate *λ*_{M} ≈ 1 day^{−1}. Using the average spheroid radii obtained from the experiments, we first obtain at steady state and then the predicted _{p} values using equation (3.8), and compare them with the marginally stable _{M,D} and _{M,S} of the shape evolution for the Darcy and Stokes models, respectively. The results are presented in tables 1⇔–3. Since the value of _{M,DS} lies between _{M,D} and _{M,S}, we do not present comparisons with the Darcy–Stokes model.

Each entry of tables 1⇔–3 shows values for the following:

— average non-dimensional radius from experiments,

— predicted

_{p}for all models,— critical

_{M,D}from the Darcy model,— critical

_{M,S}from the Stokes model, and— qualitative measure of stability from experiments.

In figures 7–9, we plot the values versus *R* for the two cases of high adhesion/low proliferation and medium adhesion/low proliferation, with the curves being the marginally stable values _{M} for the two models and the symbols corresponding to _{P}.

In table 1, we first consider the case of *l* = 3. For high adhesion and low proliferation levels, the average non-dimensional steady tumour radius is *R* = 7.5280. The predicted _{P} (square) for this case, shown in figure 7, lies above the marginally stable _{M,D} of the Darcy model and below that of _{M,S} of the Stokes model. Thus, the Stokes model predicts stable morphology, while the Darcy model predicts instability. The qualitative assessment for this case is 2 (meaning mostly stable), indicating that the Stokes prediction is more consistent with the experiments.

For medium adhesion and low proliferation levels, we have *R* = 6.656. The predicted _{P} (star) now lies above the marginal curves for both models, implying that both models predict unstable morphology, which is consistent with the qualitative measurement of a 4, corresponding to an unstable tumour. In a similar manner, we evaluate the consistency between model predictions of shape stability and experimental assessments for the rest of the table.

We conclude that the Darcy model always predicts unstable growth for the *l* = 3 mode. On the other hand, predictions by the Stokes model are consistent with the experimental data for all cases except for medium adhesion and medium proliferation. Here, the Stokes model predicts marginally unstable growth, yet the experimental assessment is a 1, implying a very stable morphology. From a mathematical standpoint, a stable tumour has a nearly spherical shape (perturbations tend to 0); the corresponding *in vitro* spheroids do not necessarily have such a shape, although they were considered most stable because of their compactness and low cell shedding. Lastly, observe that all the predicted _{P} values cluster around the _{M,S} curve for the Stokes model, suggesting that this model predicts that the experimental growth is marginally stable.

In table 2, we present comparative values for the case *l* = 4. We plot results for high adhesion with low proliferation and medium adhesion with low proliferation in figure 8. As seen in §3.1, morphologies described by the Darcy model are more sensitive to changes in *l*. Thus, the _{M,D} curve for the Darcy model shifts up considerably compared with that of the Stokes model as *l* is increased from 3 to 4. Nevertheless, the Darcy model still predicts unstable morphologies for the given experimental data, while the Stokes model remains the more predictive of the two models. In addition, owing to the insignificant shift of the _{M,S} curve, the Stokes model again predicts marginally stable growth.

Finally, we present comparative values for the *l* = 5 case in table 3. We plot results for high adhesion with low proliferation and medium adhesion with low proliferation in figure 9. Here, the upward shift of _{M,D} makes the Darcy model's prediction consistent with the experimental assessment for the high adhesion/low proliferation case, but not for the medium adhesion/low proliferation case. Now the upward shift of the curve for the Darcy model brings the marginally stable _{M,D} even closer to the predicted _{p}. The Darcy model is more predictive here than in the previous two cases of *l* = 3 and 4. Still, the Stokes curve does not shift much, and hence it predicts marginally stable tumour morphologies. Thus, the Stokes model remains more consistent with the experimental data.

## 4. Discussion

Using mathematical models of tumour growth that represent the tissue as an incompressible fluid, we studied governing factors of tumour morphological stability by assessing parameters that may control the morphology. We evaluated three single-phase sharp interface models via the following constitutive relations: Darcy's law, Stokes' law, and the combined Darcy–Stokes law. The Darcy model was previously considered by Greenspan (1976), Byrne & Chaplain (1995), Friedman & Reitich (1999), Cristini *et al.* (2003), and others, while the Stokes law was studied by Friedman & Hu (2007). Both models were considered by King & Franks (2004, 2007*a*,*b*) and Franks & King (2003). The combined Darcy–Stokes law was considered but not analysed in Zheng *et al.* (2005). Employing the Darcy model, Cristini *et al.* (2003), found that shape morphology is controlled by the competition between the ratio of apoptosis and proliferation, and the ratio of proliferation and the relaxation mechanisms (cell mobility and cell–cell adhesion). Frieboes *et al.* (2006) supported this hypothesis through experiments with *in vitro* glioblastoma cultures and obtained an experimental range of measurements of cell–cell adhesion and cell proliferation for which tumour morphology is marginally stable. Simulations of vascularized tumours *in vivo* (Cristini *et al.* 2005) showed that the morphological stability could be predicted as a function of cell substrates (e.g. oxygen and nutrients) diffusing from the vasculature, which modulate the cell adhesion and proliferation parameters. These results were successfully extended to three dimensions (using a multi-phase model, Wise *et al.* (2008)) and applied to clinical data in Bearer *et al.* (2009), providing a basis for predicting the underlying cellular conditions from the morphology observed in histopathology (e.g. hypoxic in the case of invading single cell strands or adequately vascularized in the case of rounded tumour fronts). Recent modelling work (Macklin *et al.* 2009) has highlighted the complex interactions between tumours and their vasculature, and their corresponding growth and morphology. Hydrostatic stress induced by a growing tumour can suppress blood flow in vessels by constricting their radii, thus leading to overall lower supply of cell substrates and heightened tumour morphological instability (Cristini *et al.* 2003, 2006).

We extended the analysis in Cristini *et al.* (2003) to the Stokes and Darcy–Stokes models, in order to assess which model's predictions of tumour shape stability were the most consistent with the experimental results of Frieboes *et al.* (2006). Two dimensionless parameters and were introduced that describe cell adhesion and cell death, respectively, although the definition of is model-dependent. The linear stability analyses reveal that these parameters are determinants of tumour shape morphology. We also established that the shape stability predicted by the Darcy model is more sensitive to the perturbation mode number compared with the Stokes model.

The combination of the parameter with the relative effect of vascularization (parameter *B*) subdivides tumour growth into three regimes: low, moderate, and high vascularization (Cristini *et al.* 2003). In the low vascularization regime, stable and unstable evolution occurs in all three models. On the other hand, stable growth is predicted by all three models in the high vascularization regime, indicating that vascular or tissue inhomogeneity is needed for unstable growth. However, in the moderate vascularization regime, the Stokes and Darcy–Stokes models yield both stable and unstable evolution, thus differing from the stable morphology predicted by the Darcy model.

We assessed the models in the avascular growth regime against the experimental results in Frieboes *et al.* (2006), where the growth of *in vitro* human glioblastoma spheroids was controlled by glucose (mainly moderating cell adhesion and metabolism) and growth serum (mainly moderating cell proliferation). Unlike Frieboes *et al.* (2006), we account for apoptosis in the estimation of experimental . For each particular case of glucose and growth serum, we compared model calculations of adhesion (_{p}) and the marginal stability level of adhesion (_{M,D} and _{M,S} for the Darcy and Stokes models), and assessed the model predictions of shape stability to determine model consistency against qualitative experimental observations. The results show that the Stokes model is the most consistent with the experimental results, predicting our experimental tumour growth to be marginally stable. This prediction is insensitive to the model parameters (*λ*, *l*, , *R*), and provides an indication that the model is robust. The Darcy model is overall less consistent with the experimental data, although it becomes more predictive as the polar wavenumber *l* of the perturbation is increased from 2 to 5.

In the future, we will investigate the effects of nonlinearity, using the diffuse domain approach from Li *et al.* (2009) and the numerical methods developed by Wise *et al.* (2008) to simulate the sharp interface models considered in this paper (the Darcy law was simulated earlier by Cristini *et al.* (2003) and Li *et al.* (2007) using a boundary-integral method). We will then consider other models of soft tissue mechanics such as elasto-viscoplastic constitutive laws (Ambrosi & Preziosi 2009) and develop a corresponding morphological linear stability theory.

## Acknowledgements

We thank the reviewers for their valuable comments. V.C. acknowledges partial funding from the Cullen Trust for Health Care, the Department of Defense, and the National Institutes of Health (NIH). J.L. acknowledges partial support from the NIH. V.C. and J.L. gratefully acknowledge partial support from the National Science Foundation (NSF) Division of Mathematical Sciences.

## Footnotes

- Received March 31, 2010.
- Accepted May 7, 2010.

- © 2010 The Royal Society