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 multiscale 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 threedimensional 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 selfconsistent modelling framework that can be extended to the multiscale 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 cellbased and continuum models (see recent reviews [1–19]). For example, models have been applied to brain cancer (e.g. [20–22]) and breast cancer (e.g. [23,24]). Furthermore, modelling has shown that tumour morphology may serve as a predictor of invasiveness [11,19,25–29]. In particular, Byrne & Chaplain [30,31] revealed that cell–cell adhesion and external nutrient concentration are key parameters controlling the stability of threedimensional multicellular spheroids. Extending this work, Cristini et al. [25] proposed that two nondimensional 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. [25] by employing three macroscale singlephase sharp interface models, and by applying linear stability analyses to further elucidate the determinants of tumour shape.
Greenspan [32] considered necrotic tumours in the avascular stage, where growth is regulated solely by nutrient in the surrounding microenvironment. Byrne & Chaplain [33] proposed a model for nonnecrotic tumours where nutrient is supplied through the surrounding vascularized environment. Following this work [33] and our previous work [25], we first consider nonnecrotic tumours in a vascularized environment, and later simplify the model to the avascular condition to compare its predictions to our in vitro experiments [27]. 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, [34]). 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 Cristini et al. [25], Greenspan [32], Byrne & Chaplain [33], Friedman & Reitich [35] and others, while Stokes' law was studied by Friedman & Hu [28]. Both models were investigated by King & Franks [36–38] and Franks & King [39]. The combined Darcy–Stokes law was simulated but not analysed in Zheng et al. [34]. Cristini et al. [25] showed that a nondimensionalization 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. [27] 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. [25] 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. [27] to determine the consistency between the model predictions of shape stability and the experimental data. Following [25] and using the theory of vector spherical harmonics [40], 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 singlephase 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. [27]. Briefly, the ACBT (grade IV human glioblastoma multiforme) cell line was cultured to grow tumour spheroids using a liquidoverlay technique [41]. 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 [42], thus increasing tumour cell motility as a response to hypoxia [43].
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 [33,35] and [25] and others, we present the governing equations for growth for a nonnecrotic 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.1where 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.2Hence Γ is given by 2.3where λ_{σ} is the uptake rate. We assume constant and uniform cell substrate concentration at the interface Σ (and outside Ω), 2.4
As in [32,33,35] and [25] 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.5where λ_{P} is the net cellproliferation 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 [25]: 2.7where α^{−1} is the cell motility. At the boundary, the pressure is assumed to satisfy the Laplace–Young boundary condition [44], 2.8where τ 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 [40] to the linear analysis of each model, we first recapitulate the linear results for the Darcy model in Cristini et al. [25]. We also use a different nondimensionalization than that in Cristini et al. [25].
2.3.1. Nondimensionalization
We define the following nondimensional parameters 2.11where 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 quasisteady condition in equation (2.16).
We introduce the modified concentration σ̄ such that 2.12where B is the dimensionless relative effect of vascularization, 2.13
Dropping the bars and primes, as in [25], we obtain the nondimensional Darcy system: 2.14 2.15 2.16with 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 nondimensional 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. [25] and Li et al. [45]. 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.21which depends on the apoptosis parameter and vasculature parameter B. The righthand 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.22where δ 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. [25] and Li et al. [45], the evolution equation for the shape perturbation δ/R of the interface is 2.23where 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 righthand 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 [37,39], 2.24where 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. Nondimensionalization
We rescale the stress tensor by a characteristic value T_{a}: 2.27
We use the same nondimensional 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 nondimensional parameters, 2.28where 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 nondimensional Stokes system, we have 2.29 2.30 2.31with boundary conditions 2.32 2.33 2.34where T = ∇v + (∇v)^{T} + ν(∇ · v)I − pI is the nondimensional stress tensor. The interface evolves according to the normal velocity as in equation (2.10).
The nondimensional 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 nondimensional 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.35where H(R) = (1/tanh R) − (1/R). The functions I_{l}(R) and j_{l}(ir) (where ) are the modified Bessel function of the first kind and the spherical Bessel function, respectively, related by j_{l} (ix) = 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 [46]. The term j″_{l} (iR)/j_{l} (iR) 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.37where 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 [34]. 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.39and the tumour interface evolution as in equation (2.10).
2.5.1. Nondimensionalization
We can nondimensionalize 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 nondimensional system: 2.40 2.41 2.42with boundary conditions 2.43 2.44 2.45
As before, T is the nondimensional 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 nondimensional Stokes system by the nondimensional Brinkman equation (2.41).
2.5.2. Linear stability analysis
We perform the linear stability analysis of the Darcy–Stokes model. By solving the nondimensional 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.46where 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 1a. 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 nonnegative. Hence, growth remains stable for = 0. As indicated by equation (2.23), promotes shape instability. Considering instead critical values such that dR/dt = 0 from equation (2.21), we plot _{M,D} versus the unperturbed tumour radius R for different modes l in figure 1b. 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/dt)(δ/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/dt)(δ/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/dt)(δ/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 4a, 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 4b. 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 (steadystate) value in equation (2.21), figure 4c 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 [25], 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 ̄ = /(1 − B) and cell adhesion as ̄=/(1 − B), which remains nonnegative. 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 dR/dt → −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. [26,29,47,48]), 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 [27], 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 [27], 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 dR/dt = 0, corresponding to the in vitro glioblastoma cultures when the tumour spheroids have reached a diffusionlimited steadystate size. That is, we obtain from equation (2.21) with B = 0: 3.1where R is the steadystate experimental radius scaled by the empirically observed oxygen diffusion length L ≈ 125 µm. From the experiments [27], 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 dR/dt near R ≈ 0, we obtain the net mitosis rate: 3.2
Then the evolution of the unperturbed dimensional tumour radius R̄ becomes (dR̄/dt̄) = λ̃_{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.3where 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.5where R = R̄/L is the nondimensional tumour radius at steadystate, and is the apoptosis parameter at steady state from equation (3.1).
In [27], 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.6where 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 [27], 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 nondimensional 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 nondimensional 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 singlephase 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 Cristini et al. [25], Greenspan [32], Byrne & Chaplain [33], Friedman & Reitich [35] and others, while the Stokes law was studied by Friedman & Hu [28]. Both models were considered by King & Franks [36–38] and Franks & King [39]. The combined Darcy–Stokes law was considered but not analysed in [34]. Employing the Darcy model, Cristini et al. [25], 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. [27] 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 [26] 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 multiphase model, [49]) and applied to clinical data in [50], 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 [51] 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 [25,26].
We extended the analysis in [25] 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. [27]. Two dimensionless parameters and were introduced that describe cell adhesion and cell death, respectively, although the definition of is modeldependent. 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 [25]. 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. [27], 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 [27], 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. [52] and the numerical methods developed by Wise et al. [49] to simulate the sharp interface models considered in this paper (the Darcy law was simulated earlier by Cristini et al. [25] and Li et al. [45] using a boundaryintegral method). We will then consider other models of soft tissue mechanics such as elastoviscoplastic constitutive laws [53] 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.
 Received March 31, 2010.
 Accepted May 7, 2010.
 This Journal is © 2010 The Royal Society