The intracellular and extracellular dynamics that govern tumor growth and invasiveness in vivo remain poorly understood. Cell genotype and phenotype, and nutrient, oxygen, and growth factor concentrations are key variables. In previous work, using a reaction-diffusion mathematical model based on variables that directly describe tumor cell cycle and biology, we formulated the hypothesis that tumor morphology is determined by the competition between heterogeneous cell proliferation caused by spatial diffusion gradients, e.g., of cell nutrients, driving shape instability and invasive tumor morphologies, and stabilizing mechanical forces, e.g., cell-to-cell and cell-to-matrix adhesion. To test this hypothesis, we here obtain variable-based statistics for input to the mathematical model from in vitro human and rat glioblastoma cultures. A linear stability analysis of the model predicts that glioma spheroid morphology is marginally stable. In agreement with this prediction, for a range of variable values, unbounded growth of the tumor mass and invasion of the environment are observed in vitro. The mechanism of invasion is recursive subspheroid component development at the tumor viable rim and separation from the parent spheroid. Results of computer simulations of the mathematical model closely resemble the morphologies and spatial arrangement of tumor cells from the in vitro model. We propose that tumor morphogenesis in vivo may be a function of marginally stable environmental conditions caused by spatial variations in cell nutrients, oxygen, and growth factors, and that controlling these conditions by decreasing spatial gradients could benefit treatment outcomes, whereas current treatment, and especially antiangiogenic therapy, may trigger spatial heterogeneity (e.g., local hypoxia), thus causing invasive instability. (Cancer Res 2006; 66(3): 1597-604)

Although tissue architecture and cellular environment are believed to be dominant determinants of tumor shape (13), the intracellular and extracellular factors that promote a particular tumor to adopt a specific morphology are not well understood in part because examination of tumors in vivo is difficult. The study of multicellular tumor spheroids in vitro as a model of three-dimensional in vivo tumor microenvironments may provide novel insights into tumor growth and invasion, although in vitro conditions offer the tumor cells unlimited space for expansion and lack of interaction with host cells (46). Variations in spheroid growth rates and extent of central necrosis have been attributed to fluctuations in oxygen and nutrient concentrations (7, 8) and tumor growth characteristics in vitro and in vivo have been extensively modeled (913). There have been reports that neither a limiting spheroid volume nor an increasing necrotic core could be observed for glioblastoma (14). This indicates that neither volume growth nor cell count of spheroids can be modeled by a simple exponential formulation and that spheroid volumes may fluctuate around the asymptotic volumes predicted by the Gompertz growth model (15). In fact, time-dependent oscillating patterns in tumor volumes superimposed on the Gompertz growth curve were observed for rat glioblastoma spheroids (16).

A spheroid can be viewed as a network of individual cells, each with its own proliferation potential. Growth of the spheroid is the outcome of the balance between individual cell proliferation and internal cohesive forces. Observed oscillatory growth suggests heterogeneity in tumor cells that are proliferating versus quiescent cells. Thus, spheroids may have the complexity of self-organized dynamic systems, which are regulated by both environmental and internal noise (17, 18).

Although it is usually vascular tumors that exhibit irregular shapes and complex morphology, there have been reports of irregular avascular tumors exhibiting complex growth characteristics (1922). In vitro tumors that deviate from spheroidal shape by expressing branched or chained structures have been observed (e.g., ref. 21). Spheroid cultures have been classified into three general shapes, “fused cell mass,” “tight aggregate,” and “loose aggregate” (23). Previous attempts to quantify tumor morphology include using the fractal dimension of the tumor periphery to characterize degree of tumor aggression focusing on the existence and stability of asymmetric solutions to a mathematical model of solid tumor growth (22). Thermodynamic principles applied to hepatocyte spheroid self-assembly indicate that minimization of free energy might drive morphology (24). Also, an expanding tumor may exert both mechanical pressure and traction on its microenvironment (25). The tumor mass may compress the bulk matrix radially outward but the tips of invasive cells pull the surrounding matrix inward. A feedback mechanism may link volumetric growth and invasive expansion (26). Mechanical forces influence tumor shape and may be related to tumor proliferation and invasive growth. Cellular traction may cause extracellular matrix alignment, in which matrix filaments form tracks that promote cellular elongation and directed migration, leading in some cases to formation of multicellular tubular structures (27) and intratumor cellular swirls. To explain the occurrence of invasive branches in brain tumors, a “homotype intrabranch” attraction and the principle of “least resistance, most permission, highest attraction” have been proposed (26). In this scenario, cells follow each other because of reduced mechanical resistance, enhanced haptotactic gradient, and increased chemical attraction, all as part of a self-organizing adaptive biosystem. Emergence of multicellular clusters in networks formed by migrating cells may occur as well, which may represent a guiding influence to the invasion dynamics (26).

Our results presented here support the hypothesis that tumor shape, although being affected by random cellular proliferation and adhesion and by complex mechanical interactions at the tumor viable rim, may result deterministically from the competition between individual cell proliferation and internal cohesive forces in the presence of microenvironmental substrate diffusion gradients (28). Spatial gradients in nutrient, oxygen, and growth factors levels are formed and affected by the three-dimensional heterogeneous arrangement of cells and extracellular matrix (2932), and also by the vasculature in vivo (33, 34). Intratumoral regions of hypoxia and acidosis may thus result in spatially heterogeneous tumor cell proliferation and migration (28). The need of cells to maximize their exposure to the medium to allow maximal substrate uptake may lead to shape instability and invasive tumor morphology [i.e., in low-substrate medium, cells strive to minimize cell contact (invasive morphology) and maximize exposure of all of the membrane transporters to media]; where substrate is abundant, cells may be better served by cell-to-cell contact (compact spheroid morphology). Under these conditions, invasive morphology manifests itself with the development of low-wave-number deterministic fluctuations in cell positions at the tumor viable rim as cells proliferate and regulate adhesion based on local diffusion gradients, leading to recursive formation and growth of a finite number of subtumors that eventually may break off from the parent tumor. The growth of low-wave-number modes is characteristic of diffusion-driven instabilities in materials and biomaterials (see refs. 35, 36 and references therein). We use an in silico model (37, 38) to quantify this instability. We show that this computational model can be accurately fitted to in vitro experimental data to determine variables, and that multidimensional computer simulations based on these input variables predict tumor cell spatial arrangement and tumor morphologies that closely resemble those observed in vitro. Tumor morphologies can thus be quantified via a predictive mathematical model with the long-term goal of optimizing variables for therapy application to minimize growth and invasion. Diffusion-driven tumor shape instability could be suppressed in vivo by enforcing more spatially homogeneous nutrient and oxygen supply because normoxic conditions act both by decreasing gradients and increasing cell adhesion and, consequently, the mechanical forces that maintain well-defined tumor morphology (28).

Cell culture. ACBT (human glioblastoma) and BT4C (rat glioma) cell lines were cultured in DMEM (Life Technologies Invitrogen, Carlsbad, CA) supplemented with 10% fetal bovine serum (FBS; Life Technologies Invitrogen), 2 mmol/L l-glutamine, and 1% penicillin/streptomycin in humidified 7.5% CO2 at 37°C. At 70% confluence, cells were removed from incubation and left at room temperature for ∼20 minutes. The resultant cell clusters (∼10 cells) were transferred to a Petri dish precoated with 1% agar in DMEM basal medium and grown to tumor spheroids using a liquid-overlay technique (39). After 1 to 2 days in culture, the spheroids were removed and placed individually in 48-well plates (Corning, Costar 3548) coated with 0.75% agar (FMC Bioproducts, Rockland, ME). Medium was replenished every 3 days.

Variation in nutrient and growth factors. For both cell lines, three sets of multiwell plates were used, one for each level of glucose. DMEM with high (4.5 g/L) and low glucose (1 g/L) were purchased as standard basal medium. Medium glucose (2.75 g/L) was prepared by mixing the other two. Each set consisted of three individual multiwell plates for 1%, 5%, and 10% FBS concentrations.

Spheroid measurements. Spheroid growth was observed using a Leitz microscope at magnification ×100 (11 observations per plate for 36 days). Spheroid diameter was measured via the arithmetic average of two orthogonal variables, one being the longest observable axis. The field of vision was round with a diameter of 1,800 μm. Photographs were taken with an Olympus camera at magnification ×160 mounted on top of the microscope with a photography window of 1,130 × 1,430 μm. Spheroids were fixed in 4% paraformaldehyde and embedded in paraffin. Histologic cross-sections were obtained by slicing in 6 μm increments and staining according to H&E.

Color-coded staining. Before generating cell clusters, one half of cell population was membrane labeled with the PKH26 red fluorescent cell linker kit and the other half with PKH67 green (both from Sigma-Aldrich). Spatial distribution of red- and green-labeled cells in spheroids was determined using a Zeiss Fluorescent microscope.

Reaction-diffusion model of tumor progression. Briefly, we model (37, 38) nutrient, oxygen, and growth factor diffusion (Eq. A) through tumor interstitium and uptake by tumor cells (40); local specific tumor mass growth rate (Eq. B) as the divergence of velocity u of tumor cells; local pressure p in tissue by Darcy's law (Eq. C1) with boundary condition (Eq. C2) at tumor-host tissue interface (41):

\[0\ =\ L^{2}{\nabla}^{2}\ \frac{\mathit{n}}{\mathit{n}_{\mathrm{a}}}\ {-}\ \frac{\mathit{n}}{\mathit{n}_{\mathrm{a}}};\]
\[\mathit{{\lambda}}_{\mathrm{M}}^{{-}1}{\nabla}\ {\bullet}\mathrm{\mathbf{u}}\ =\ \frac{\mathit{n}}{\mathit{n}_{\mathrm{a}}}\ {-}\ A;\]
\[\mathrm{\mathbf{u}}\ =\ {-}\mathit{{\lambda}}_{\mathrm{M}}\mathit{L}^{2}{\nabla}\frac{\mathit{p}}{\mathit{p}_{\mathrm{a}}};\]
\[\frac{\mathit{p}}{\mathit{p}_{\mathrm{a}}}\ =\ G^{{-}1}\mathit{L{\kappa}}\]

where n is local nutrient, oxygen or growth factor concentration, whichever dominates growth; na is medium concentration; pa = λML2 / μ is characteristic pressure magnitude in the tissue; diffusion length L = (D / η)1/2, where η is the rate of nutrient uptake within tumor cells (inverse time) and D is diffusion coefficient; and λM is rate of accumulation of tumor cell mass due to cell mitosis (inverse time). Tumor cells and extracellular matrix are treated as comprising porous media with hydraulic conductivity μ that is a measure of cell mobility; A = λD / λM is a dimensionless variable, where λD is a death rate (inverse time) that describes disintegration of cell mass and radial effusion of intracellular fluid away from necrotic regions (41), thus controlling the overall tumor size (∇ • u = 0 outside the tumor tissue). For simplicity, the “tumor” phase encompasses tumor cell matter, interstitial fluid, and extracellular matrix. Cell adhesion forces throughout the spheroid are modeled using an equivalent surface tension γ at the tumor/tissue interface (of local curvature κ). The dimensionless variable G = λM / λR controls morphological stability, where λR = μγL−3 is a relaxation rate (inverse time) associated with cell adhesion. For simplicity, we also assume that viable cell mass density is uniform in the tumor (41) and that regions become necrotic where nutrient and oxygen concentration decrease below some specified minimum (38). The above reaction-diffusion model is solved using an adaptive (42, 43) finite-element/level-set method (38) in two spatial dimensions.

The three-dimensional linear stability analysis presented in ref. (37) for the reaction-diffusion mathematical model briefly described in Materials and Methods revealed that mass growth and shape evolution of in vitro spheroids depend on four microphysical variables: cell mitosis rate λM (inverse time), diffusion length L, and two dynamics-controlling dimensionless variables G (ratio of tumor cell mitosis rate to mechanical relaxation rate, e.g., associated to cell adhesion) and A (ratio of cell death rate to mitosis rate). In the model formulations, cell mitosis and death rates describe production and destruction of tumor cell mass associated with these processes. The hypothesis underlying this theory is that spheroid morphology is controlled by the variable G describing the competition between shape-destabilizing forces linked to proliferation and stabilizing, relaxation forces linked primarily to cell-to-cell and cell-to-matrix adhesion. For G small, spheroids are morphologically stable (i.e., spherical or nearly spherical) and mass growth and steady spheroid radius following a period of growth are functions of the variable A only, regulating the balance between mass production through mitosis and mass destruction through cell death. For the spheroids in our cell culture, the equations governing the growth of the spheroid equivalent radius R (radius of a sphere of equal volume, rescaled with L) and the spheroid morphologic stability (t is time rescaled with λM−1) are as follows:

Mass growth:

\[\frac{\mathrm{d}\mathit{R}}{\mathrm{d}\mathit{t}}\ =\ {-}A\ {\cdot}\ R/3\ +\ 1/\mathrm{tanh}\ R\ {-}\ R^{{-}1}\]

Morphologic stability5

Note that Eq. E is incorrect in ref. (37).

:

\[A\ =\ 54G^{{-}1}\ {\cdot}\ R^{{-}3}\ {-}\ \frac{\mathrm{3}}{\mathrm{4}}\ +\ \frac{\mathrm{3}}{\mathrm{4}}\left(\mathit{7R}^{{-}\mathrm{1}}\ +\ \frac{\mathit{I}_{\mathrm{11/2}}(\mathit{R})}{\mathit{I}_{\mathrm{9/2}}(\mathit{R})}\right)(\mathrm{1/tanh}\ \mathit{R\ {-}\ R}^{{-}\mathrm{1}})\]

(the I's are modified Bessel functions) as can be determined from a three-dimensional linear stability analysis (37) of the model described in Materials and Methods.

We calculate the microphysical variables by fitting the above model to in vitro data of spheroid growth. These variables can then provide a measure of tumor growth and morphologic stability. Experimental observations (see Materials and Methods) consisted of nine cases: three levels of glucose, each with three different concentrations of FBS. Typical spheroid cross-sections and curves of spheroid diameter versus time are reported in Fig. 1. Spheroids showed a rim of viable, proliferating cells surrounding a hypoxic core with predominantly necrotic cells. As expected (4446), spheroid growth curves showed initial exponential growth [early growth rate from Eq. D is dR / d

\(\overline{t}\)
λMR / 3, where λM
\(\overline{t}\)
t, thus the initial slope in the log-linear plot in Fig. 1 (left) is λM / 3] and eventual plateau reaching a stationary radius (dR / dt = 0). Diffusion gradients established along the radial distance from the spheroid center limit cell viability to a rim across which levels of oxygen and glucose drop below the cell viability limit. Thus, from measurement of the thickness of the viable rim (Fig. 1, right), we estimate a diffusion length

Figure 1.

Left, growth history of glioma spheroids in medium glucose level (see Materials and Methods) at different FBS concentrations. Dashed line, closed points, 1% FBS; dashed line, open points, 5% FBS; solid line, closed points, 10% FBS. (Corresponding data for different levels of glucose are qualitatively similar although the mitosis rate, λM, increases slightly with glucose level.) Right, typical spheroid histologic cross-section showing proliferating cells at the viable rim surrounding a hypoxic core. Bar, 130 μm.

Figure 1.

Left, growth history of glioma spheroids in medium glucose level (see Materials and Methods) at different FBS concentrations. Dashed line, closed points, 1% FBS; dashed line, open points, 5% FBS; solid line, closed points, 10% FBS. (Corresponding data for different levels of glucose are qualitatively similar although the mitosis rate, λM, increases slightly with glucose level.) Right, typical spheroid histologic cross-section showing proliferating cells at the viable rim surrounding a hypoxic core. Bar, 130 μm.

Close modal
\[\mathit{L}\ {\approx}\ 100\ \mathrm{{\mu}}\mathit{m}.\]

By fitting Eq. D to the growth curves in Fig. 1, we obtain an average proliferation rate

\[\mathit{{\lambda}}_{\mathrm{M}}\ {\approx}\ 1\mathrm{day}^{{-}1}\]

and a range of values for the death variable:

\[\mathrm{0.26}\ {\leq}\ \mathit{A}\ {\leq}\ \mathrm{0.38.}.\]

Note that Eq. D predicts that final spheroid sizes (obtained by setting the time derivative dR / dt = 0) are dependent only on the ratio A of death-to-mitosis rates. This is consistent with our in vitro observations where the terminal size of all spheroids is roughly the same regardless of glucose and serum levels (see, e.g., Fig. 1, left). Note also that the death rate, assumed to be uniform to simplify the mathematical analysis, may be spatially heterogeneous. For example, it will presumably be highest at the perinecrotic region between the hypoxic core and the viable rim. Computer simulations shown below correspond to a more sophisticated model accounting for spatially nonuniform death by necrosis (A nonuniform).

According to this model, during spheroid growth, cell mitosis occurs in the viable rim where nutrients, oxygen, and growth factor levels are adequate, and cell death (predominantly by necrosis) occurs in the inner regions where diffusion limitations prevent these substances from being present in adequate levels. Figure 1 (left) indicates that during the first 3 weeks, tumors grow in size. Correspondingly, there is an unbalance of mass fluxes: in the hypoxic core, dead cells disintegrate releasing intracellular fluid that is pushed outwards and contributes to form new cells in the proliferating rim; at the interface of proliferating and hypoxic regions, cells die and are pushed into the necrotic core by the cells inside the rim that require more volume as they are proliferating, thus creating an inward flux of dead cells. After ∼1 month, the accumulation of dead cells in the hypoxic core balances cell proliferation in the viable rim and, thus, these two opposite mass fluxes equilibrate and the overall spheroid growth is arrested. It is instructive to directly estimate the value of the equilibrium variable A from a scaling argument based on these considerations. The differential equation (Eq. B) in Materials and Methods describes local conservation of mass in the spheroid (37). By using the steady profile of cell nutrients or oxygen, n / na = (R / sinh R) r−1 sinh r (r is radial spherical coordinate measured from the center of the spheroid and all lengths have been rescaled with the diffusion length L), obtained by solving the reaction-diffusion equation (Eq. A) in the spheroid (37) and integrating Eq. B over the entire spheroid volume gives the rate of change of the total spheroid mass (Eq. D). By setting dR / dt = 0 for a stationary spheroid, and solving for A, we obtain A = 3R−2 (R cosh R − sinh RR) / sinh R. Note that n / na = 1 at the spheroid surface (r = R) and decreases inside the spheroid and approaching the center of the spheroid. From inspection of the stationary spheroids, we have a typical dimensionless stationary radius R ≈ 8 (i.e., a diameter of 1,600 μm). This and the above formula give the estimate A ≈ 0.33, consistent with the results of direct fits to the growth curves.

A rigorous quantitative analysis of how our experimental conditions translate into different values of the microphysical model variable G controlling the morphologic stability of spheroid growth would require more sophisticated, targeted experiments to directly measure how cell adhesion and proliferation depend on growth factor, glucose, and oxygen levels and is beyond the scope of the present investigation. Instead, here, to estimate the variable G, we develop the above scalings further. The differential equation (Eq. C1; Darcy's law) describes the formation of local pressure gradients in the tissue (due to proliferation) and the following cell motion down these gradients, thus promoting mass growth and host invasion by the tumor cells. The boundary condition (Eq. C2; Laplace-Young law) describes the balance of such pressure and adhesion forces, where pa is a characteristic value of pressure inside the proliferating rim. From Eq. C2, by setting ppa / 3 (i.e., ∣u∣ ≈ M / 3) and estimating, for a morphologically stable spheroidal tumor of radius

\(\overline{\mathit{R}}\)
(dimensionless radius R =
\(\overline{\mathit{R}}\)
/L), the surface curvature κ ≈ 2
\(\overline{\mathit{R}}\)
−1, we predict: G ≈ 6L/
\(\overline{\mathit{R}}\)
= 6R−1. From the statistics of stationary spheroid radius and thickness of viable rim in our in vitro experiments, we obtain the range of values where morphologically stable spheroids can exist:

\[0:6\ {\leq}\ \mathit{G}\ {\leq}\ 0:9\]

Spheroids with values of G above this range will not be stable as cell adhesion forces are too weak. Note that the above scaling provides an indirect method for estimating values of G necessary for stability using the measurements of spheroid radii and without the need of directly measuring the strength of cell adhesion. This is because this quantity does not affect the final size of morphologically stable spheroids. We will later show that the above condition is necessary but is not sufficient to guarantee morphologic stability. Note also that in Eq. C1, the surface tension force that here describes cell adhesion is proportional to the local surface curvature. Thus, according to this model, adhesion forces will strongly prevent fingering morphologies (corresponding to highest curvature values) but less strongly oppose the development of low-wave-number instabilities (characterized by low curvatures).

Figure 2 is a morphologic stability diagram obtained by plotting the curve of stationary spheroid radius (from Eq. D with dR / dt = 0), and the curves from Eq. E controlling morphologic stability, over the ranges of values for the variables A and G estimated from in vitro statistics (shaded area) of morphologically stable spheroids. The part of the curve for Eq. D that lies within the experimental range of values of A denotes possible stationary spheroid radii R. For a given G, the region of the plane on the left of the corresponding curve from Eq. E describes morphologically stable spheroids, whereas the region on the right describes unstable spheroids (37). Higher values of G correspond to curves for Eq. E shifted to the left, thus reducing the occurrence of morphologically stable spheroids. We conclude that for given values of A and G, a spheroid of stationary radius R will remain morphologically stable only if its position on the curve for Eq. D is on the left of the stability curve for Eq. E corresponding to its value of G. Accordingly, the spheroid denoted by the filled square is very stable, being on the left of all curves for Eq. E 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). Because our in vitro spheroids exhibited stationary radii R ≈ 8, Fig. 2 predicts that most of these tumor spheroids are only marginally stable and shape instability may develop. Note that in vitro (and in vivo), different regions of a tumor may have different values of G and A and, thus, some regions may be stable, whereas others may be unstable.

Figure 2.

Morphologic stability diagram showing death variable A versus spheroid radius R (rescaled with diffusion length L). The curves for given values of G are obtained from plotting Eq. E governing spheroid morphologic stability. Experimental conditions for morphologically stable spheroids (shaded area) are enclosed by these and by the horizontal lines delimiting the range of values of A, all estimated by fitting the mathematical model to the in vitro data. The curve “stationary radius” is obtained by setting dR / dt = 0 in Eq. D governing spheroid radius growth. Three representative stationary spheroid radii are reported on this curve as sampled in vitro. Because this curve crosses and continues beyond the shaded region, we conclude that most glioma spheroids under these in vitro conditions are marginally stable.

Figure 2.

Morphologic stability diagram showing death variable A versus spheroid radius R (rescaled with diffusion length L). The curves for given values of G are obtained from plotting Eq. E governing spheroid morphologic stability. Experimental conditions for morphologically stable spheroids (shaded area) are enclosed by these and by the horizontal lines delimiting the range of values of A, all estimated by fitting the mathematical model to the in vitro data. The curve “stationary radius” is obtained by setting dR / dt = 0 in Eq. D governing spheroid radius growth. Three representative stationary spheroid radii are reported on this curve as sampled in vitro. Because this curve crosses and continues beyond the shaded region, we conclude that most glioma spheroids under these in vitro conditions are marginally stable.

Close modal

In our in vitro experiments, we varied glucose and serum concentrations to vary cell adhesion and rates of cell proliferation. The corresponding range of values of the variable G was estimated above for morphologically stable spheroids. Rates of cell proliferation increased with serum concentration as expected. Cases with 1% FBS had slowest tumor mass growth, whereas 10% FBS cases had the fastest (see Fig. 1, left, and Eq. B in Materials and Methods with n being the local concentration of growth factors in serum). 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 reduce oxygenation in tumor spheroids >1 mm in diameter (47) and that hypoxia can increase tumor cell motility by increased production of autocrine motility factor and by up-regulation of hepatocyte growth factor (4851). Thus, an important effect of higher glucose was to decrease cell adhesion forces and thereby contribute to morphologic instability. Most stable, compact morphologies (0.6 ≤ G ≤ 0.9) were observed with low/medium levels of both glucose and FBS, in which tumors 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 (G > 0.9). 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 (G > 0.9). According to our theory, the heterogeneous access to oxygen, nutrients, and growth factors at the cellular scale due to spatial gradients of these substances established by diffusion led to differential proliferation within the viable rim of the in vitro tumor spheroids, with groups of cells proliferating faster in regions of higher substance concentrations. Under these conditions, fluctuations in shape were created as cells proliferated heterogeneously at the viable rim and would lead, for high G (i.e., relatively weak cell adhesion), to morphologic instability manifesting itself with the development and growth of subspheroidal structures. Sometimes, these subspheroids, composed of groups of predominantly viable cells, would break off completely from the “mother” spheroid and grow into separate tumor spheroids (see photographs in Figs. 3 and 4). Thus, unstable glioma spheroids grew as a composite of spheroidal subtumors. This process repeated itself on the subtumor, leading to recursive subspheroid growth as the main mechanism of tumor spheroid morphogenesis. Unstable morphologies in vitro seemed to exhibit mostly tumor surface perturbations characterized by low wave numbers (e.g., 3 or 4) at the onset of the instability, consistent with our diffusional instability theory (37). Accordingly, Eq. E was adapted from the theory of ref. (37) for a wave number equal to 4. This roughly means that three or four “bumps” are expected to develop on the spheroid surface and eventually lead to subspheroids (see photographs in Figs. 3 and 4). In contrast, for low G, surface perturbations would not grow and the overall spheroidal tumor shape would remain stable. In conclusion, we were able to observe in the experiments both stable and unstable spheroids by varying the variable G according to our theory.

Figure 3.

Spheroid morphologies from computer simulations and experiments. Low-wave-number instabilities arise on spheroid surfaces eventually leading to the development and separation of subspheroids. Simulation snapshots (length rescaled with diffusion length L, time t rescaled with mitosis time λM−1 ≈ 1 day) showing the outer boundary and inner perinecrotic rims (A), and local levels of diffusing substances (B), such as oxygen or glucose. Photographs C, glioma spheroids growing in culture. Subspheroids are highlighted in middle and bottom photographs. Bar, 130 μm.

Figure 3.

Spheroid morphologies from computer simulations and experiments. Low-wave-number instabilities arise on spheroid surfaces eventually leading to the development and separation of subspheroids. Simulation snapshots (length rescaled with diffusion length L, time t rescaled with mitosis time λM−1 ≈ 1 day) showing the outer boundary and inner perinecrotic rims (A), and local levels of diffusing substances (B), such as oxygen or glucose. Photographs C, glioma spheroids growing in culture. Subspheroids are highlighted in middle and bottom photographs. Bar, 130 μm.

Close modal
Figure 4.

Left, cell displacement velocities (arrows) during subspheroid formation (zoom from simulation of Fig. 3 at t = 30). Right (top and bottom), histologic cross-sections showing cell arrangements in subspheroids grown in vitro. In the photographs, subspheroids can be seen forming from parent spheroids. As a result of shape instability on the parent spheroid, swirling patterns of cell arrangements can be identified surrounding necrotic cells in the subspheroids both in the simulations and the experiments. Bar, 130 μm.

Figure 4.

Left, cell displacement velocities (arrows) during subspheroid formation (zoom from simulation of Fig. 3 at t = 30). Right (top and bottom), histologic cross-sections showing cell arrangements in subspheroids grown in vitro. In the photographs, subspheroids can be seen forming from parent spheroids. As a result of shape instability on the parent spheroid, swirling patterns of cell arrangements can be identified surrounding necrotic cells in the subspheroids both in the simulations and the experiments. Bar, 130 μm.

Close modal

We did computer simulations of glioma spheroid growth in two spatial dimensions using the numerical algorithm developed in ref. (38) to solve in the nonlinear regime (large shape deformations) the mathematical model briefly summarized in Materials and Methods. Compact spheroid morphology (data not shown) was achieved for values of G within the stable region (as defined in Fig. 2). For an unstable case (large G >0.9), snapshots of the evolution of a tumor spheroid are shown in Fig. 3 as predicted from a computer simulation. The outer boundary tracks the surface of the spheroid; the inner boundary encloses regions of hypoxia, where necrotic cells are also present. A thin rim (thickness roughly equal to L) of viable and actively proliferating cells is predicted, surrounding a large hypoxic core (compared with Fig. 1, right). Note the irregularities arising on the spheroid surface from random oscillations in the positions of the cells. These oscillations introduce low-wave-number perturbations on the spheroid surface, as can be seen in the snapshots from the simulation and in the photographs from the in vitro experiments. These perturbations grow (large G), leading to formation of subspheroidal structures that eventually separate from the mother spheroid (see for comparison the photographs from the in vitro experiments in Fig. 3). Clusters of spheroids are thus formed that allow the tumor mass to grow to a much larger size and over a much larger region than would have been possible had the spheroid maintained a compact shape and instability had not occurred, as predicted in ref. (28). If the spheroid remains compact, nutrient and oxygen diffusion limitations to mass growth cause the spheroid to reach a final, stable size (see Fig. 1). In the last snapshot from our simulation in Fig. 3, it is shown that the instability repeats itself on the subspheroids. This was also observed in the experiments.

Our computational results and experimental in vitro observations of morphology are in agreement, thus supporting the theory that invasive tumor morphologies observed in vitro (and universally in vivo) are the result of a shape instability driven by diffusion gradients in the tumor microenvironment. This morphologic instability develops as a result of differential cell proliferation induced by underlying diffusion gradients across the tumor mass where cells that are exposed to higher levels of nutrients and oxygen are favored for proliferation. Perturbations of cell positions on the surface of a tumor or of spatial distribution of these substances trigger the instability. This differential growth becomes more pronounced as time proceeds because cells at the leading edge of a shape perturbation (“bulb” of cells, e.g., see Fig. 3) have higher exposure to nutrient, thereby gaining a proliferative advantage, whereas cells at the receding edge become increasingly disadvantaged. In Fig. 3, the gradients of concentration n / na of the substance that is most important to proliferation (either oxygen, glucose, or growth factors) are shown as predicted from the simulation. These gradients drive the instability, and proliferation of the subspheroids follows the gradients. We conclude that morphologic instability provides a very effective means to invasion because it allows tumor cells better access to oxygen and nutrients in the surrounding environment (28).

In Fig. 4, a detail of subspheroid formation is shown from the simulation and from the experiment of Fig. 3. In the simulation snapshot, mass fluxes (associated to the velocity u described in Materials and Methods) are depicted by arrows. Outside the tumor, these are associated to the flow of fluid in the culture dish driven by proliferation and motion of cells. In the tumor, the arrows represent cell mass fluxes. As cells proliferate, they are pushed by the rising pressure p as described by Eq. C1. Cell mass in the unstable regions (proliferating rim, acting as a source of cell mass) protrudes based on outward movement of cells at the leading edge and undergoes involution in the area of “pinching” distal to this growth (hypoxic regions, acting as a sink of cell mass), leading to eventual tumor breakup and separation of subspheroids from the parent spheroid. The resulting swirling patterns are typical of motion in the presence of a source and a sink in fluids and materials. Swirling arrangements of viable cells are also roughly identifiable in the photographs, as is commonly observed experimentally with fibroblast-like cells. Note that the flux of cells into the hypoxic regions predicted by the simulation is exaggerated in magnitude as an artifact of the model used (described in Materials and Methods) that considers a uniform cell mass density. In reality, cell mass density in the hypoxic regions is much lower than in the viable rim because of necrosis; hence, the real flux of cells is also lower.

The aggressive proliferation and invasiveness of tumor metastases (and often of primary tumors) lead to high levels of patient mortality. An improved understanding of morphologic patterns during tumor growth may aid design of therapies aimed at limiting or neutralizing their invasiveness and recurrence after therapeutic intervention. With this goal in mind, we investigated the growth and shape of human gliomas in vitro and in silico as a function of oxygen, nutrient and growth factor concentrations and of cell adhesion forces. To test the hypothesis that invasive infiltrative morphology may be a function of diffusion gradients (28), we varied glucose and serum concentrations for human and rat glioblastoma spheroids in vitro. Parameters obtained from these experimental observations were then used in a multiscale and multidimensional linear analysis and computer simulation of tumor growth (37, 38) to predict, and possibly explain, glioma spheroid morphology and invasiveness. Results presented here indicate that diffusion gradients of nutrients and growth factors may have a direct bearing on morphology and invasion. Tumor shape is a function of cell proliferation and adhesion based on cellular access to oxygen, nutrient and growth factors. When proliferation rates are relatively low, cell adhesion is sufficient to maintain compact tumor shapes. When proliferation rates are higher, subtumors may emanate from the main tumor, separating off from the tumor surface as clusters of cells. This situation maximizes the overall cell population's access to nutrient by exposing more tumor surface area to the outside environment. (Note that in vivo this effect may be mitigated as it also exposes more tumor cell surface area and associated antigens to immune attack.) Recursive subtumor growth can be mathematically represented by a morphologic instability triggered by low-wave-number fluctuations at the tumor viable periphery. The underlying causes for this instability are diffusion gradients of nutrients and growth factors in the tumor microenvironment. This mechanism may help to explain the morphology of an aggressive genotype recently observed for breast cancer in vitro with transformed mammary epithelial cells (52), where excess proliferation and luminal filling produce a complicated structure characterized by multiple acinar structures (53).

One of the main assumptions in this study is that tumor spheroid growth in vitro and in silico represents, to some extent, tumor growth in vivo, and that tumor morphogenesis, as a function of access to nutrient and growth factors, occurs similarly in vitro and in vivo. The results presented here reflect avascular tumor growth in the absence of effects from host cells, which could modify these results. Tumors in vivo possess varying degrees of vascularity and have a complex interface with normal cells, including immune cells that may be toxic, and the extracellular matrix, which may be quite different from that generated in vitro. Also, tumors in vivo are profoundly acidotic with pH typically below 7 and often below 6.4. A full model of tumor invasiveness must incorporate all of these variables. We have described tumor invasion as the result of a morphologic instability. An alternative explanation for our in vitro observations may be that recursive subtumor growth is mainly a function of individual cell proliferation, such as stem cells, as was recently suggested (54). Although this conjecture deserves further investigation, our observations imply that under the experimental conditions in this study, tumor growth was promoted equally at the tumor surface by all cells, which combined to form subtumor structures. To check this, we did experiments in which half of the cells in a spheroid were stained in one color and the other half in another color (see Materials and Methods). The subspheroids formed afterward were observed to have cells in both colors, thus suggesting that they did not originate from one individual cell with higher proliferative potential than its neighbors.

By using a sophisticated computer simulation of tumor growth, we established a mathematical framework that enabled prediction and interpretation of experimental results, showing that unstable tumor morphologies are driven by nutrient spatial variations. Our results support the hypothesis (28) that heterogeneous oxygen and nutrient supply may drive tumor growth and invasiveness through a diffusional instability mechanism. These considerations may be relevant during chemotherapy, radiotherapy, and antiangiogenic therapy, all of which could introduce spatial (and temporal) variations in oxygen and cell nutrients (55). In contrast, tumors may approach a compact noninvasive morphology when cell adhesion or other stabilizing mechanical forces (e.g., tumor encapsulation) are maximized. We propose that compact tumor morphologies may be achievable by maintaining uniform nutrient levels at the cellular scale and homogeneous microenvironmental conditions, thus suppressing instability (28).

We are currently working to test this hypothesis in an in vivo tumor xenograft model. Human breast and prostate cancer growth is observed using a window chamber following s.c. implantation in a mouse model. As the implanted tumor grows and is vascularized, the increase and growth of long-wave perturbations at the tumor/host interface is observed, as in the in vitro and in silico results presented here. This unstable morphology well correlates with measurements of steep radial gradients of pH and oxygen pointing away from the tumor.

Note: This article is dedicated to the memory of Rutherford Aris.

Grant support: NIH Laser Microbeam and Medical Program (LAMMP, P41-RR01192; B. Tromberg and C-H. Sun), NIH grants U56CA113004 and R01 CA093650 from the National Cancer Institute (R. Gatenby), and National Science Foundation and National Cancer Institute (V. Cristini and H.B. Frieboes).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

We thank Steen J. Madsen (University of Nevada, Las Vegas, NV) for help in planning and setting up the experiments, and Henry Hirschberg (Rikshospitalet, Oslo, Norway), John Lowengrub, and Xiangrong Li (University of California Irvine) for helpful comments and insights.

1
Sansone BC, Delsanto PP, Magnano M, Scalerandi M. Effects of anatomical constraints on tumor growth.
Phys Rev E
2001
;
64
:
021903
.
2
Schmeichel KL, Weaver VM, Bissell MJ. Structural cues from the tissue microenvironment are essential determinants of the human mammary epithelial cell phenotype.
J Mammary Gland Biol Neoplasia
1998
;
3
:
201
–13.
3
Enam SA, Rosenblum ML, Edvardsen K. Role of extracellular matrix in tumor invasion: migration of glioma cells along fibronectin-positive mesenchymal cell processes.
Neurosurgery
1998
;
42
:
599
–608.
4
Sutherland RM. Cell and environment interactions in tumor microregions: the multicell spheroid model.
Science
1988
;
240
:
177
–84.
5
Hamilton G. Multicellular spheroids as an in vitro tumor model.
Cancer Lett
1998
;
131
:
29
–34.
6
Santini MT, Rainaldi G, Indovina PL. Apoptosis, cell adhesion and the extracellular matrix in three-dimensional growth of multicellular tumor spheroids.
Crit Rev Oncol Hematol
2000
;
36
:
75
–87.
7
Freyer JP, Sutherland RM. Determination of diffusion constants for metabolites in multicell tumor spheroids.
Adv Exp Med Biol
1983
;
159
:
463
–75.
8
Freyer JP, Sutherland RM. Regulation of growth saturation and development of necrosis in EMT6/Ro multicellular spheroids by the glucose and oxygen supply.
Cancer Res
1986
;
46
:
3504
–12.
9
Demicheli R, Pratesi G, Foroni R. The exponential-Gompertzian growth model: data from six tumor cell lines in vitro and in vivo. Estimate of the transition point from exponential to Gompertzian growth and potential clinical applications.
Tumori
1991
;
77
:
189
–95.
10
Rygaard K, Spang-Thomsen M. Quantitation and Gompertzian analysis of tumor growth.
Breast Cancer Res Treat
1997
;
46
:
303
–12.
11
Ferrante L, Bompadre S, Possati L, Leone L. Parameter estimation in a Gompertzian stochastic model for tumor growth.
Biometrics
2000
;
56
:
1076
–81.
12
Adam JA. Mathematical models of prevascular spheroid development and catastrophe—theoretic description of rapid metastatic growth/tumor remission.
Invasion Metastasis
1996
;
16
:
247
–67.
13
Marusic M, Bajzer Z, Freyer JP, Vuk-Pavlovic S. Analysis of growth of multicellular tumour spheroids by mathematical models.
Cell Prolif
1994
;
27
:
73
–94.
14
Nirmala C, Rao JS, Ruifrok AC, Langford LA, Obeyesekere M. Growth characteristics of glioblastoma spheroids.
Int J Oncol
2001
;
19
:
1109
–15.
15
Chignola R, Schenetti A, Andrighetto G, et al. Forecasting the growth of multicell tumour spheroids: implications for the dynamic growth of solid tumours.
Cell Prolif
2000
;
33
:
219
–29.
16
Chignola R, Schenetti A, Chiesa E, et al. Oscillating growth patterns of multicellular tumour spheroids.
Cell Prolif
1999
;
32
:
39
–48.
17
Coffey DS. Self-organization, complexity and chaos: the new biology for medicine.
Nat Med
1998
;
4
:
882
.
18
Wiesenfeld K, Moss F. Stochastic resonance and the benefits of noise: from ice ages to crayfish and SQUIDS.
Nature
1995
;
373
:
33
.
19
Bredel-Geissler A, Karbach U, Walenta S, Vollrath L, Mueller-Klieser W. Proliferation-associated oxygen consumption and morphology of tumor cells in monolayer and spheroid culture.
J Cell Physiol
1992
;
153
:
44
–52.
20
Mueller-Klieser W. Three-dimensional cell cultures: from molecular mechanisms to clinical applications.
Am J Physiol Cell Physiol
1997
;
273
:
C1109
–23.
21
Enmon RM, Jr., O'Connor KC, Lacks DJ, Schwartz DK, Dotson RS. Dynamics of spheroid self-assembly in liquid-overlay culture of DU 145 human prostate cancer cells.
Biotechnol Bioeng
2001
;
72
:
579
–91.
22
Byrne H, Matthews P. Asymmetric growth of models of avascular solid tumours: exploiting symmetries.
IMA J Math Appl Med Biol
2002
;
19
:
1
–29.
23
Hedlund TE, Duke RC, Miller GJ. Three-dimensional spheroid cultures of human prostate cancer cell lines.
Prostate
1999
;
41
:
154
–65.
24
Tzanakakis ES, Hansen LK, Hu WS. The role of actin filaments and microtubules in hepatocyte spheroid self-assembly.
Cell Motil Cytoskeleton
2001
;
48
:
175
–89.
25
Gordon VD, Valentine MT, Gardel ML, et al. Measuring the mechanical stress induced by an expanding multicellular tumor system: a case study.
Exp Cell Res
2003
;
289
:
58
–66.
26
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 spheroid model.
Cell Prolif
2001
;
34
:
115
–34.
27
Vernon RB, Angello JC, Iruela-Arispe ML, Lane TF, Sage EH. Reorganization of basement membrane matrices by cellular traction promotes the formation of cellular networks in vitro.
Lab Invest
1992
;
66
:
536
–47.
28
Cristini V, Frieboes HB, Gatenby R, Caserta S, Ferrari M, Sinek J. Morphological instability and cancer invasion.
Clin Cancer Res
2005
;
11
:
6772
–9.
29
Acker H, Carlsson J, Mueller-Klieser W, Sutherland RM. Comparative pO2 measurements in cell spheroids cultured with different techniques.
Br J Cancer
1987
;
56
:
325
–7.
30
Carlsson J, Acker H. Relations between pH, oxygen partial pressure and growth in cultured cell spheroids.
Int J Cancer
1988
;
42
:
715
–20.
31
Casciari JJ, Sotirchos SV, Sutherland RM. Glucose diffusivity in multicellular tumor spheroids.
Cancer Res
1988
;
48
:
3905
–9.
32
Filion RJ, Popel AS. A reaction-diffusion model of basic fibroblast growth factor interactions with cell surface receptors.
Ann Biomed Eng
2004
;
32
:
645
–63.
33
Jain RK. Determinants of tumor blood flow: a review.
Cancer Res
1988
;
48
:
2641
–58.
34
Jain RK. Physiological barriers to delivery of monoclonal antibodies and other macromolecules in tumors.
Cancer Res
1990
;
50
:
814
–9s.
35
Cristini V, Lowengrub J. Three-dimensional crystal growth-I: linear analysis and self-similar evolution.
J Cryst Growth
2002
;
240
:
267
–76.
36
Cristini V, Lowengrub J. Three-dimensional crystal growth-II: nonlinear simulation and suppression of the Mullins-Sekerka instability.
J Cryst Growth
2004
;
266
:
552
–67.
37
Cristini V, Lowengrub J, Nie Q. Nonlinear simulation of tumor growth.
J Math Biol
2003
;
46
:
191
–224.
38
Zheng X, Wise S, Cristini V. Nonlinear simulation of tumor necrosis, neo-vascularization and tissue invasion via an adaptive finite-element/level-set method.
Bull Math Biol
2005
;
67
:
211
–59.
39
Sutherland RM, Carlsson J, Durand RE, Yuhas J. Spheroids in cancer research.
Cancer Res
1981
;
41
:
2980
–94.
40
Byrne H, Chaplain M. Growth of nonnecrotic tumors in the presence and absence of inhibitors.
Math Biosci
1995
;
130
:
151
–81.
41
Chaplain MAJ. Avascular Growth, Angiogenesis and vascular growth in solid tumours: the mathematical modelling of the stages of tumour development.
Math Comput Model
1996
;
23
:
47
–87.
42
Cristini V, Blawzdziewicz J, Loewenberg M. An adaptive mesh algorithm for evolving surfaces: simulations of drop breakup and coalescence.
J Computational Phys
2001
;
168
:
445
–63.
43
Anderson A, Zheng X, Cristini V. Adaptive unstructured volume remeshing—I: The method.
J Comput Phys
2005
;
208
:
616
–25.
44
Gompertz B. On the nature of the function expressive of the law of human mortality, and on a new mode of determining the value of life contingencies.
Philos Trans R Soc Lond
1825
;
115
:
513
–85.
45
Winsor CP. The Gompertz curve as a growth curve.
Proc Natl Acad Sci U S A
1932
;
18
:
1
–8.
46
Laird AK. Dynamics of tumor growth.
Br J Cancer
1964
;
18
:
490
–502.
47
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.
48
Young SD, Marshall RS, Hill RP. Hypoxia induces DNA overreplication and enhances metastatic potential of murine tumor cells.
Proc Natl Acad Sci U S A
1988
;
85
:
9533
–7.
49
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.
50
Cairns RA, Kalliomaki T, Hill RP. Acute (cyclic) hypoxia enhances spontaneous metastasis of KHT murine tumors.
Cancer Res
2001
;
61
:
8903
–8.
51
Postovit LM, Adams MA, Lash GE, Heaton JP, Graham CH. Oxygen-mediated regulation of tumor cell invasiveness. Involvement of a nitric oxide signaling pathway.
J Biol Chem
2002
;
277
:
35730
–7.
52
Debnath J, Brugge J. Modelling glandular epithelial cancers in three-dimensional cultures.
Nat Rev Cancer
2005
;
5
:
675
–88.
53
Debnath J, Mills K, Collins N, Reginato M, Muthuswamy S, Brugge J. The role of apoptosis in creating and maintaining luminal space within normal and oncogene-expressing mammary acini.
Cell
2002
;
111
:
29
–40.
54
Singh SK, Clarke ID, Terasaki M, et al. Identification of a cancer stem cell in human brain tumors.
Cancer Res
2003
;
63
:
5821
–8.
55
Sinek J, Frieboes H, Zheng X, Cristini V. Two-dimensional chemotherapy simulations demonstrate fundamental transport and tumor response limitations involving nanoparticles.
Biomed Microdev
2004
;
6
:
297
–309.