Abstract
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 reactiondiffusion 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., celltocell and celltomatrix adhesion. To test this hypothesis, we here obtain variablebased 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): 1597604)
 computer simulation
 tumor invasion
 tumor spheroids
 tumor morphology instability
 mathematical model
Introduction
Although tissue architecture and cellular environment are believed to be dominant determinants of tumor shape ( 1– 3), 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 threedimensional 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 ( 4– 6). 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 ( 9– 13). 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, timedependent 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 selforganized 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 ( 19– 22). 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 selfassembly 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 selforganizing 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 threedimensional heterogeneous arrangement of cells and extracellular matrix ( 29– 32), 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 lowsubstrate 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 celltocell contact (compact spheroid morphology). Under these conditions, invasive morphology manifests itself with the development of lowwavenumber 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 lowwavenumber modes is characteristic of diffusiondriven 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 longterm goal of optimizing variables for therapy application to minimize growth and invasion. Diffusiondriven 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 welldefined tumor morphology ( 28).
Materials and Methods
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 lglutamine, and 1% penicillin/streptomycin in humidified 7.5% CO_{2} 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 liquidoverlay technique ( 39). After 1 to 2 days in culture, the spheroids were removed and placed individually in 48well 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 crosssections were obtained by slicing in 6 μm increments and staining according to H&E.
Colorcoded 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 SigmaAldrich). Spatial distribution of red and greenlabeled cells in spheroids was determined using a Zeiss Fluorescent microscope.
Reactiondiffusion 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 tumorhost tissue interface ( 41):where n is local nutrient, oxygen or growth factor concentration, whichever dominates growth; n_{a} is medium concentration; p_{a} = λ_{M}L^{2} / μ 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 reactiondiffusion model is solved using an adaptive ( 42, 43) finiteelement/levelset method ( 38) in two spatial dimensions.
Results
The threedimensional linear stability analysis presented in ref. ( 37) for the reactiondiffusion 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 dynamicscontrolling 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 shapedestabilizing forces linked to proliferation and stabilizing, relaxation forces linked primarily to celltocell and celltomatrix 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:
Morphologic stability 5:(the I's are modified Bessel functions) as can be determined from a threedimensional 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 crosssections 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 ( 44– 46), spheroid growth curves showed initial exponential growth [early growth rate from Eq. D is dR / d ≈ λ_{M}R / 3, where λ_{M} ≈ t, thus the initial slope in the loglinear 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 By fitting Eq. D to the growth curves in Fig. 1, we obtain an average proliferation rate and a range of values for the death variable: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 deathtomitosis 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 / n_{a} = (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 reactiondiffusion 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 R − R) / sinh R. Note that n / n_{a} = 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; LaplaceYoung law) describes the balance of such pressure and adhesion forces, where p_{a} is a characteristic value of pressure inside the proliferating rim. From Eq. C2, by setting p ≈ p_{a} / 3 (i.e., ∣u∣ ≈ Lλ_{M} / 3) and estimating, for a morphologically stable spheroidal tumor of radius (dimensionless radius R = /L), the surface curvature κ ≈ 2 ^{−1}, we predict: G ≈ 6L/ = 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: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 lowwavenumber 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.
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 upregulation of hepatocyte growth factor ( 48– 51). 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.
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 lowwavenumber 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 / n_{a} 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 fibroblastlike 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.
Discussion
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 lowwavenumber 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 longwave 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.
Acknowledgments
Grant support: NIH Laser Microbeam and Medical Program (LAMMP, P41RR01192; B. Tromberg and CH. 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.
Footnotes

Note: This article is dedicated to the memory of Rutherford Aris.
 Received September 2, 2005.
 Revision received November 14, 2005.
 Accepted November 28, 2005.
 ©2006 American Association for Cancer Research.