PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of plosonePLoS OneView this ArticleSubmit to PLoSGet E-mail AlertsContact UsPublic Library of Science (PLoS)
 
PLoS One. 2010; 5(11): e13942.
Published online 2010 November 30. doi:  10.1371/journal.pone.0013942
PMCID: PMC2994713

Emergent Properties of Tumor Microenvironment in a Real-Life Model of Multicell Tumor Spheroids

Terence Lee, Editor

Abstract

Multicellular tumor spheroids are an important in vitro model of the pre-vascular phase of solid tumors, for sizes well below the diagnostic limit: therefore a biophysical model of spheroids has the ability to shed light on the internal workings and organization of tumors at a critical phase of their development. To this end, we have developed a computer program that integrates the behavior of individual cells and their interactions with other cells and the surrounding environment. It is based on a quantitative description of metabolism, growth, proliferation and death of single tumor cells, and on equations that model biochemical and mechanical cell-cell and cell-environment interactions. The program reproduces existing experimental data on spheroids, and yields unique views of their microenvironment. Simulations show complex internal flows and motions of nutrients, metabolites and cells, that are otherwise unobservable with current experimental techniques, and give novel clues on tumor development and strong hints for future therapies.

Introduction

Multicellular tumor spheroids (MTS) stand out as the most important in vitro model of pre-vascular solid tumors [1][8]. MTS often have a regular, almost spherical structure, and their apparent simplicity has led to repeated attempts to capture their features with neat mathematical models. However, the absence of vascularization and the near sphericity hide an internal complexity which is not easy to tame either with analytic mathematical models [9][12], or with numerical models based on rough simplifications of the biological settings such as cellular automata or other lattice-based models [13][16]. Moreover the presence of a growing necrotic core [1] and of an extracellular matrix [17], the appearance of convective cell motions [18], and the heterogeneous response to chemotherapics [19], point to the importance of MTS as an in vitro model of tumors, and most of all to their relevance to understand tumor heterogeneity, but they also point to the difficulties of producing a useful, predictive model of MTS.

The appearance of widely different resistance phenomena to antitumor therapies in similarly grown, isolated MTS of the same cell type [19] indicates that random fluctuation phenomena play an all-important role in the growth kinetics of MTS. It is well-known that the discrete events at the single-cell level (like transitions from one cell-cycle phase to the next, mitosis, cell death, etc.) do display some randomness, and one can pinpoint the source of large-scale variability on these fluctuations, as they are amplified and propagated by cell-cell and cell-environment interactions. Thus, the complexity of MTS development can only be captured by a fine-grained, multiscale model, and we need a mathematical description at the single-cell level. Since cells communicate with other cells and the environment, the other actors of this complex play are the concentration gradients of important molecular species that depend on the structure of the extracellular space and of the facilitated transport processes into and out of individual cells, and the mechanical forces that push and pull cells as they proliferate with repeated mitoses and then shrink after death [20]. These processes mix with complex nonlinear interactions between the biochemical and the mechanical part, and this highlights again the importance of an effective model at the single-cell level.

On the basis of such motivations, we have developed a numerical model of MTS that incorporates a working model of single cells [21], [22]. We have first put forward a broad outline of its structure in reference [23], and it differs from other models developed in the past [9][16] because it captures at the same time both the basic features of cell metabolism, growth, proliferation and death, and provides a true lattice-free calculation of cell motions, as they are pushed and pulled by the forces exerted by dividing cells, by the growth of other cells, and by the shrinking of dead cells. We also wish to stress that the model parameters are either derived from experiment or are deduced from reasonable theoretical arguments, so that, essentially, there are no free parameters – there can only be some residual variability in biophysically meaningful ranges – the model is truly predictive, and the results are not merely qualitative but quantitative as well.

Here we illustrate in broad terms the structure of the program and report the results of the first simulations of single spheroids (technical implementation details are relegated to Text S1). We find that the simulations agree quite well with experimental measurements on real spheroids, and show unexpected and important internal patterns. Moreover, we wish to stress that the methods delineated in this paper represent very general practical solutions to problems that are common to any simulation of cell clusters, and they are just as important.

Biochemical behavior of individual cells

The elementary building blocks in this model of MTS are the individual tumor cells that behave as partly stochastic automata [21], [22]. Figure 1 summarizes the biochemical pathways that are included in the single-cell model: cell metabolism is driven by oxygen, glucose and glutamine, and transforms these substances into energy molecules, molecular building blocks and waste products, following the well-known biochemical reaction chains [24]. Further details can be found in the original papers [21], [22] and in Text S1, which also includes important upgrades to the original model [21], [22].

Figure 1
Rough sketch of the biochemical pathways incorporated in the model of single cells.

In the present version of the program, the stochasticity is mostly concentrated in the discrete events: for instance, mitochondria are partitioned at random between daughter cells at mitosis, and cells can die because of metabolite accretion, according to a Poissonian cytotoxicity model (see Text S1).

We remark that in this approach glutamine also stands for the wider class of aminoacids, and lactate is the paradigm of all metabolites: we use the concentrations of glutamine and lactate to represent these two classes of substances in phenomenological parameterizations wherever needed. Similarly we use the number of mitochondria and ATP content to model the dynamics of cell volume; the single-cell model also contains representative members of the cyclin protein class to compute the passage of checkpoints and entry into the different cell phases [21], [22], [25], [26], and finally into mitosis (see also figure S1 for a sketch of the cell cycle in the simulation program).

The complete map of the biochemical pathways included in the simulation program is shown in figure S2. This map comprises only the most basic pathways, however we cannot afford to introduce a more complex network at this stage of program development. Indeed, our final aim is the simulation of MTS with a volume as large as 1 mmAn external file that holds a picture, illustration, etc.
Object name is pone.0013942.e001.jpg, which corresponds to more than one million cells, so that simulation results overlap actual experimental measurements [19], [27], [28]. Since the differential system involves 19 independent biochemical variables per cell, we must eventually integrate at least 19 million coupled nonlinear differential equations for the biochemical cell variables alone (this grows to at least 25 million equations when we include the position and velocity variables), and thus even this minimal single-cell model leads to a daunting computational task (see Text S1 for further details on the algorithmic complexity of the program).

Reaction-diffusion processes and the environment

Substances like oxygen are transported into and out of cells by normal diffusion while molecules like glucose require facilitated diffusion processes. This means that cell membranes play an important role for substances like glucose, and that in this case the diffusion of each such molecular species towards cells in the inside of a spheroid needs the free volume in the extracellular space to proceed, and that we must model this space as well as the cells to obtain a realistic simulation. We have shown how to do this in reference [29], where we have also discussed ways to tame the exceedingly high stiffness of the very large set of reaction-diffusion and transport equations that arise in this context (see also Text S1). The external environment itself is included in these equations, and evolves in time as well. In the present model, each cell contributes 15 internal variables and 4 extracellular variables: these extracellular variables are the masses of oxygen, glucose, glutamine and lactate in the extracellular volume surrounding the cell. Because of its smallness, the extracellular space has an extremely short characteristic filling time, which can be as fast as few tens of microseconds. On the other hand, the macroscopic features of MTS evolve over times as long as months (i.e., times of the order of An external file that holds a picture, illustration, etc.
Object name is pone.0013942.e002.jpg), and thus the numerical integrator must be able to handle phenomena that span 12 orders of magnitude in time [29]. The internal biochemical reactions included in the numerical model are much slower and their fastest characteristic times are only as low as An external file that holds a picture, illustration, etc.
Object name is pone.0013942.e003.jpg, much longer than the diffusion times [29], [30]. The topology of diffusion in the extracellular spaces is obviously dictated by the cells themselves, and the program uses the network of cells centers as the scaffolding for the corresponding discretized diffusion problem. The links between the cells' centers – i.e., the proximity relations – are provided by a Delaunay triangulation [31], [32], which is computed repeatedly [33] as the cluster of cells grows and rearranges itself under the pushes and pulls of volume growth, mitosis, and the shrinking of dead cells (see also figure S3). Moreover, the proliferation of cells means that both the number of cells and the total number of links steadily grow, and that the differential system of equations that model metabolism, transport and diffusion changes all the time, and becomes increasingly complex. The 3D Delaunay triangulation itself is not an exceedingly heavy computational burden for the program, as it turns out that efficient algorithms can compute it, on average, with An external file that holds a picture, illustration, etc.
Object name is pone.0013942.e004.jpg time computational complexity [33][35], so that this algorithm is indeed feasible for very large clusters of cells.

Biomechanical evolution of the simulated MTS

Real cells have passive viscoelastic mechanical features, but they also move actively under the pushes of their own cytoskeleton, and to the best of our knowledge there is no comprehensive model of cellular biomechanics [36], [37]. Thus, we resort once again to phenomenological simplifications, and the first and foremost is that our cells are stretchable spheres, characterized by their radius, and a few other parameters that specify their viscoelastic properties (see Text S1 for a more detailed description and the list of parameters). We also specify a pairwise interaction force between cells, repulsive when a cell pushes against a neighbor, and attractive when we try to detach it from its neighbor. For small deviations from the equilibrium distance, we assume that the interaction force is described by the Hertz model (explained in Text S1), while for large deformations due to compression we set the force to a fixed saturation value, and for large distances the attractive force decays to zero (see figure S4). The description of the interaction forces is tuned to hold also during mitosis (see Text S1 and figure S5). Even though this is a rough approximation of the overall mechanical behavior of cells, there are many details that must be managed to make it work, and they are all described in Text S1.

Here the Delaunay triangulation that we use as the scaffolding for the diffusion problem turns out to be useful once again: the same cell-cell links also define the set of neighbors of each cell, and therefore the global problem of computing the pairwise interactions between cells can be reduced to a single loop over all cells and the small limited number of their immediate neighbors, so that this operation has an An external file that holds a picture, illustration, etc.
Object name is pone.0013942.e005.jpg computational complexity only – and it does not grow when we include the cost of the Delaunay triangulation [35] – instead of the An external file that holds a picture, illustration, etc.
Object name is pone.0013942.e006.jpg complexity of generic pairwise interactions.

Results

The first and most obvious result is the outstanding match of the growth curves of simulated spheroids with those of real spheroids: figure 2 shows a few stages of the growth of a simulated spheroid (a real spheroid is shown for comparison in figure 3), while figure 4 compares the growth curve of a single simulated spheroid with the growth curves of real spheroids grown in vitro. Here we see that the growth curves are very much alike, and we found that simulation runs with different parameters – in the biophysically meaningful ranges – produce very similar growth curves, in spite of structural internal changes: the growth curves are thus rather robust with respect to parameter changes.

Figure 2
Snapshots of one simulated spheroid taken at different times.
Figure 3
Photograph of a spheroid grown in vitro from HeLa cells in agar.
Figure 4
Growth curve of a simulated tumor spheroid (solid line).

Several experiments [37][42] have yielded many accurate measurements of oxygen and glucose concentrations and other quantities vs. spheroid radius; these values are part of the output of our simulation program as well (see figure 5 and figure 6), and a comparison with the experimental data is shown in table 1. On the whole the agreement of simulation data of single spheroids with the experimental values is quite good, and we wish to stress that this is not the result of a fit a posteriori, but rather of the a priori choice of model features and parameters. These results qualify as true predictions of the numerical model.

Figure 5
Concentrations in the simulated spheroid.
Figure 6
Plots of the normalized average intracellular concentration of lactate (green), glucose (blue), and ATP (red).
Table 1
Comparisons with experimental parameters.

The necrotic core of spheroids is another important feature that is well reproduced in the simulations, and it is clearly visible in central slices of the simulated spheroid in figure 2. The simulations also provide detailed, quantitative snapshots of the necrotic core dynamics; the left column of figure 7 shows the percentage of dead cells vs. distance from the centroid of a simulated spheroid at different times. In these snapshots we can clearly observe the formation of the sharp step that marks the edge of the necrotic core.

Figure 7
Fraction of dead cells (left column) and average radial velocity (right column) at different times.

These results indicate that the simulation program is reliable and robust and reproduces – both quantitatively and qualitatively – known experimental results. However, it yields much more than just successful comparisons: figure 8 shows two views of the spheroid microenvironment that at present would be unobtainable by other means at this level of resolution. The left panel of figure 8 is a plot of the flow of glucose in the extracellular spaces of a mature spheroid, superposed on a density plot of extracellular glucose concentration, and it shows – rather unexpectedly – that there is an outward flow of extracellular glucose from the central necrotic region. In the external, viable rim the flow is inward bound, and there is a spherical shell where the flow is stationary. The right panel of figure 8 shows the corresponding plot of cell velocities in the same central slice, and we see that the velocity vectors point outward in the viable rim, while there are well-formed vortices in the central region, and the region in-between displays distinctive chaotic motions: these three regions closely match the three regions in the left panel. The right column in figure 7 shows radial velocity vs. distance from the centroid of the simulated spheroid, and sheds some more light on the nature of this structure: as more and more cells die and the necrotic core forms, the dead cells shrink and the core contracts. The contraction of the necrotic core expels the residual glucose in the extracellular spaces and produces the observed outward flow. We found that this behavior is strongly dependent on the particular value of the diffusion coefficient and on the metabolic activity of cells. In some simulations – where we used a lower value for the effective diffusion coefficient of oxygen – we observed a similar structure with oxygen as well. We remark that in the case of lactate we found no such structure, and we obtained a pH value – derived from the distribution of lactate inside the spheroid – that is very close to experimental measurements: this indicates that the discretized reaction-diffusion scheme used in the simulation program performs correctly, and that the observed flows are not algorithmic artifacts.

Figure 8
Two views of the microstructure of a simulated spheroid, with about An external file that holds a picture, illustration, etc.
Object name is pone.0013942.e033.jpg diameter and 296264 cells (183893 live cells+112371 dead cells).

Discussion

Although the program described in this paper is based on a model of individual cells that includes only the basic cell functions, the simulation results compare very well with experimental measurements, and give strong hints on the sources of individual spheroid variability. Moreover, the images obtained in single runs reveal unexpected and interesting correlations and an elaborate structure of the tumor microenvironment that could never be observed before. This unexpected, complex microstructure – the formation of different regions, and the flows that characterize them, along with the complex velocity field – can be discerned in the flows of the other substances, though not all of them, according to their effective diffusion coefficient and their metabolism: the figures of these flows are shown at full-size as supporting information. Thus if we suppose that, in a more complete description, there are An external file that holds a picture, illustration, etc.
Object name is pone.0013942.e034.jpg substances that characterize the spheroid microenvironment, and assume that the spherical shell that divides the two main regions lies in the same position for all of these substances and that their effective diffusion coefficients are uncorrelated, then An external file that holds a picture, illustration, etc.
Object name is pone.0013942.e035.jpg different spheroid structures are determined by diffusion alone. The variation of some critical parameter (e.g., a slight change in the metabolic activity due to local fluctuations in the number of dead cells, and thus a change in the effective diffusion coefficients) can potentially act as a switch and determine widely different fates for similar spheroids. This variability cannot be discerned from growth experiments: the simulations that we have performed to date indicate that the growth curve alone is not enough to distinguish between such different states, because it does not change much even when important substances, like oxygen, diffuse in markedly different ways. These different states represent different biochemical configurations of tumor microenvironment, that might exert distinct selective pressures on cells during tumor evolution.

The spheroid microstructure that is well evidenced in figure 8, and in figures S8, S9, S10, S11, S12, S13, S14, S15, S16, S17, S18, S19, S20, S21 and in Movies S1, S2, S3, shows highly correlated fluctuations that produce, e.g., islets of proliferating cells in the sea of dead cells of the core, and cell and mass flows that follow preferential channels. There is a sort of spheroid-specific self-organization of the internal structure due to these correlated fluctuations. Similar cell flows have been observed in the lab and a recent review has stressed the great significance of such findings [43]: the simulations suggest that the whole topic of cell flows and extracellular diffusion should be investigated further. On the basis of the simulation results, we also conjecture that the flow of therapeutic drugs may be diverted as well, and let some viable, proliferating tumor cells escape treatment. This means that the simulation program could eventually become an important tool to design novel treatment schedules, and possibly validate the effects of anti-tumor drugs.

Certainly the model is far from complete, and we plan to add soon several new features, like a basic model of intracellular acidity, now accounted for by a simple phenomenological parameterization, and the effects of pH and salt concentration on diffusion. However, already in its present form, we believe that this numerical model is a true testbed of biological complexity and a real virtual laboratory, and also a source of important biomedical clues.

Methods

The simulation program is written in ANSI C++: this programming language was a natural choice from the very start for distinct reasons:

  • C++ is an object-oriented language, and in a simulation such as this, it is very natural to define objects that have a clear-cut biological meaning;
  • at present, C++ programming is supported by a vast array of scientific libraries, and this helps reducing program development time;
  • the availability of the flexible and powerful C++ library CGAL [33] that handles the computational geometry structures utilized by the program (convex hulls, Delaunay triangulations and alpha shapes);
  • the availability of powerful development tools and highly optimized compilers.

The structure of the program reflects the organization explained in the paper: a layout is shown in figure 9. The functional blocks work as follows:

Figure 9
Functional blocks of the simulation program.

Initialization

At start, the environmental concentrations are set at their standard levels (see Text S1), and internal variables of all cells are set at approximate standard values (see Text S1 for the cells' variables and the physical values that are hard-coded in the program). During initialization, cells are allowed to grow and proliferate freely in an environment that is held fixed. The number of cells is also kept constant, and when a mitosis occurs one of the daughter cells is discarded. In this initial phase cells can have large oscillations of their metabolic parameters, and can occasionally step in parameter regions that would normally spell death: this does not occur here. Initialization lasts until the oscillations of metabolic parameters die out. We have determined the duration of the initialization phase observing the desynchronization of a population of initially synchronized cells: when oscillations of the relative fractions of cells in each cell-cycle phase become undetectable we estimate that cells have reached a stable state. It turns out that a simulated time of An external file that holds a picture, illustration, etc.
Object name is pone.0013942.e036.jpg (i.e. about 35 days of simulated time) is sufficient for initialization of cell with a period of about 20 hours. Usually the starting number of cells is quite small (normally just one cell to seed the growth of a single spheroid), and initialization executes in very short real time (a few seconds).

Metabolism, diffusion, transport, and growth

This part of the program solves the combined differential system of equations that describe internal cell metabolism and diffusion in the extracellular spaces (described in detail in Text S1), using the implicit Euler method. This leads to a system of nonlinear equations, that are solved in turn with a variant of the Newton-Raphson method. The functional scheme of this important part of the program is shown in figure 10. We wish to stress that although the number of variables can be quite large (more than An external file that holds a picture, illustration, etc.
Object name is pone.0013942.e037.jpg loop variables), convergence is reasonably fast, because the initial concentration values are invariably very close to the final ones.

Figure 10
Functional blocks of the C++ method that computes metabolic and extracellular variables.

Cell motion

Cell motion is also regulated by differential equations and the solution uses a strategy based on a semi-implicit method (described in detail in Text S1). Volume growth is regulated by the part that handles metabolism and diffusion, therefore it is loosely coupled to cell motion. However we have implemented an updating mechanism that effectively decouples the two parts of the program: this means that the program can use multithreading with shared memory and exploit the features of multicore processors.

Cellular events

This part of the program handles discrete events, like cell-cycle transitions, mitosis and cell death. In case of mitosis it also initializes the daughter cells – using the metabolic variables of the mother cell – and allocates memory for the new cells.

Geometry and topology of cell cluster

Geometrical and topological informations are updated here, using calls to CGAL methods [33] that compute the convex hull of the cluster of cells, the Delaunay triangulation of cell centers, and the alpha shape of the cluster – with an alpha parameter [33] equal to An external file that holds a picture, illustration, etc.
Object name is pone.0013942.e039.jpg where An external file that holds a picture, illustration, etc.
Object name is pone.0013942.e040.jpg is the average cell radius. This part of the program uses this basic information to set all relevant geometrical and topological variables in the program.

Summary statistics and dump on file

The last step in the loop computes several statistics and outputs them on summary files. It also writes periodically the whole configuration of cells on file for further processing.

Program termination

The program repeats the loop until one of the stop conditions is met: either all cells are dead, or the program executed the required number of steps. Text S1 contains additional considerations on algorithmic complexity and on measured performance (see also figure S6 and figure S7).

Additional processing to extract useful informations from the simulation data is performed with several standard tools, like Mathematica [44].

Supporting Information

Text S1

Including tables and additional references.

(0.89 MB DOC)

Figure S1

Sketch of the cell phases included in the simulation program. The arrow lengths suggest the relative duration of each phase. Phase G1 is divided in two parts: an initial subphase where cells are sensitive to variations in the environmental nutrient concentration (G1m) and a later subphase where cells are insensitive to deprivation of nutrients (G1p); in between these subphases there is an energetic checkpoint [1][3].

(0.19 MB TIF)

Figure S2

Sketch of the metabolic network. Variables within circles represent molecular species and are expressed in units of concentration or mass. Non-obvious symbols are as follows (the suffixes ext and int denote, respectively, extracellular and intracellular variables): G = glucose, G6P = glucose-6-phosphate, STORE = glucose stored in the form of glycogen, AL = lactic acid, A = glutamine, ATPp = pool of ATP molecules, DNA = nuclear mass of DNA (normalized to 1 for the whole genome), mtDNA = mitochondrial DNA. Rates are represented by squares. The red circuit represents the oxygen sensor, whereas the green circuit represents the ATP sensor [1], [2]. “Cell cycle checkpoints” denotes the molecular circuit of cell cycle control that has been modeled on the basis of previous studies on the dynamics of the allosteric effect [19], [20]. The biological foundations of this simplified metabolic network have been given in references 1 and 2. Recent improvements with respect to our previous model include: internalization rates of glucose, glutamine and lactate are sensitive to extracellular pH and this dependence is now described by smoothed functions (see the text for details); synthesis of cellular proteins, nuclear DNA and mitochondrial DNA are now described by double-substrate Michaelis-Menten chemical reactions to take into account the dependence of protein and nucleic acid biosynthesis on glutamine (which stands phenomenologically for the wider class of aminoacids) and ATP availability.

(0.77 MB TIF)

Figure S3

The geometry and topology of diffusion. a) Most substances are carried into and out of the cell by facilitated diffusion and there is an active mass exchange between cell and extracellular space. Each cell in the simulation program has its own extracellular space. b) Extracellular spaces are interconnected and there is a diffusion flow through the network of connections. c) The network of interconnected spaces is defined by a Delaunay triangulation. In this 2D representation, for any red dot we can define a Voronoy cell, i.e., the set of points in the plane that are closer to the given dot than to any other dot in the set. The dual structure is the Delaunay triangulation (Voronoy cells are black and Delaunay links are green). There is a Delaunay link between any two dots only if the respective Voronoy cells touch each other, therefore we can use the Delaunay triangulation to define proximity. This enables us to set up a discretized version of diffusion between extracellular spaces, like in b). In addition to the topology of contacts between cells we also keep into account geometry: gab in part b) is a geometric factor that modulates diffusion. d) The actual simulation is in 3D: here the Delaunay triangulation of a small cluster of cell centers shows up in transparency.

(1.37 MB TIF)

Figure S4

Pictorial representation of the interaction force between two cells. The solid curve shows qualitatively the behavior of the interaction force, while the insects depict the corresponding situations (A. cells are compressed against each other, force is repulsive; B. cells are in contact, the total force vanishes; C. cell centers are slightly apart, force is attractive because of adhesive molecules on the cells' membranes; D. cells are no longer in contact, the total force vanishes again). The inset on the right corner shows the definitions of the basic geometric variables.

(0.36 MB TIF)

Figure S5

The geometry of mitosis. R0 is the radius of the initial cells, while R1 and R2 are the radii of the daughter cells: because of random asymmetries during mitosis, the daughter cells usually have different sizes. The program places the two daughter cells inside the region initially occupied by the mother, and the axis connecting the centers points in a random direction. This forces the two cells to push one against the other, and as cells separate the axis rotates, so that the new configuration fits the positions of neighboring cells. The cells' centers are separated by a distance which is roughly 0.4 R0.

(0.16 MB TIF)

Figure S6

CPU time needed to simulate 1 hour, vs. the number of cells in the spheroid. In this run, the precision of the solution of the global diffusion transport and metabolism problem is fixed at 1% and the timestep is 50 s (so that tCPU is actually the CPU time needed to simulate 72 elementary timesteps). The solid curve shows the fit (S.69) in Text S1.

(0.15 MB TIF)

Figure S7

Total CPU time vs. N. This figure shows the total CPU time vs. N in the same run as figure S6. The solid curve is a simple fit with a quadratic polynomial function.

(0.14 MB TIF)

Figure S8

Oxygen concentration and flow at simulated time = 14 days. Cells are coloured to show the oxygen concentration (color mapping, blue = low concentration, red = high concentration), while the yellow arrows show the oxygen flow (arrow length proportional to flow intensity).

(2.06 MB TIF)

Figure S9

Oxygen concentration and flow at simulated time = 16 days. Cells are coloured to show the oxygen concentration (color mapping, blue = low concentration, red = high concentration), while the yellow arrows show the oxygen flow (arrow length proportional to flow intensity).

(1.14 MB TIF)

Figure S10

Oxygen concentration and flow at simulated time = 17 days. Cells are coloured to show the oxygen concentration (color mapping, blue = low concentration, red = high concentration), while the yellow arrows show the oxygen flow (arrow length proportional to flow intensity).

(4.88 MB TIF)

Figure S11

Extracellular glucose concentration and flow at simulated time = 14 days. Cells are coloured to show the concentration of extracellular glucose (color mapping, blue = low concentration, red = high concentration), while the yellow arrows show the flow of extracellular glucose (arrow length proportional to flow intensity).

(2.05 MB TIF)

Figure S12

Extracellular glucose concentration and flow at simulated time = 15 days. Cells are coloured to show the concentration of extracellular glucose (color mapping, blue = low concentration, red = high concentration), while the yellow arrows show the flow of extracellular glucose (arrow length proportional to flow intensity).

(2.88 MB TIF)

Figure S13

Extracellular glucose concentration and flow at simulated time = 16 days. Cells are coloured to show the concentration of extracellular glucose (color mapping, blue = low concentration, red = high concentration), while the yellow arrows show the flow of extracellular glucose (arrow length proportional to flow intensity).

(1.15 MB TIF)

Figure S14

Extracellular glucose concentration and flow at simulated time = 17 days. Cells are coloured to show the concentration of extracellular glucose (color mapping, blue = low concentration, red = high concentration), while the yellow arrows show the flow of extracellular glucose (arrow length proportional to flow intensity).

(4.97 MB TIF)

Figure S15

Extracellular glutamine concentration and flow at simulated time = 14 days. Cells are coloured to show the concentration of extracellular glutamine (color mapping, blue = low concentration, red = high concentration), while the yellow arrows show the flow of extracellular glutamine (arrow length proportional to flow intensity).

(2.06 MB TIF)

Figure S16

Extracellular glutamine concentration and flow at simulated time = 17 days. Cells are coloured to show the concentration of extracellular glutamine (color mapping, blue = low concentration, red = high concentration), while the yellow arrows show the flow of extracellular glutamine (arrow length proportional to flow intensity).

(4.94 MB TIF)

Figure S17

Extracellular glutamine concentration and flow at simulated time = 18 days. Cells are coloured to show the concentration of extracellular glutamine (color mapping, blue = low concentration, red = high concentration), while the yellow arrows show the flow of extracellular glutamine (arrow length proportional to flow intensity).

(6.11 MB TIF)

Figure S18

Extracellular lactate concentration and flow at simulated time = 18 days. Cells are coloured to show the concentration of extracellular lactate (color mapping, blue = low concentration, red = high concentration), while the yellow arrows show the flow of extracellular lactate (arrow length proportional to flow intensity). Lactate always flows outward in the simulation: this is a single snapshot taken after both glucose and glutamine have developed their split flow regime.

(6.07 MB TIF)

Figure S19

Velocity vectors projected on the plane of the slice at simulated time = 14 days. The vectors show the cells' motions in the plane of the slice (yellow arrows, arrow length proportional to flow intensity). Cells in the core perform complex looping motions, while cells in the viable rim always push outward.

(1.62 MB TIF)

Figure S20

Velocity vectors projected on the plane of the slice at simulated time = 16 days. The vectors show the cells' motions in the plane of the slice (yellow arrows, arrow length proportional to flow intensity). Cells in the core perform complex looping motions, while cells in the viable rim always push outward.

(3.08 MB TIF)

Figure S21

Velocity vectors projected on the plane of the slice at simulated time = 18 days. The vectors show the cells' motions in the plane of the slice (yellow arrows, arrow length proportional to flow intensity). Cells in the core perform complex looping motions, while cells in the viable rim always push outward.

(5.00 MB TIF)

Movie S1

Movie of a central slice of a simulated tumor spheroid showing the development of the necrotic core (red = live cells, black = dead cells).

(5.99 MB MOV)

Movie S2

Movie of a central slice of a simulated tumor spheroid showing the flow of extracellular glucose (same coding as figures S11, S12, S13, S14).

(21.18 MB MOV)

Movie S3

Movie of a central slice of a simulated tumor spheroid showing the map of projected cell velocities (same coding as figures S19, S20, S21).

(28.40 MB MOV)

Footnotes

Competing Interests: The authors have declared that no competing interests exist.

Funding: This work has been supported by the Italian Institute for Nuclear Physics (INFN) under the names VIRTUS, VBL and VBL-Rad. Relevant information can be found on the INFN website, http://www.infn.it/indexen.php. The funder had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

References

1. Sutherland RM. Cell and environment interactions in tumor microregions: the multicell spheroid model. Science. 1988;240:177–84. [PubMed]
2. Bjerkvig R. Spheroid Culture in Cancer Research. Boca Raton, Fla.: CRC Press; 1992.
3. Mueller-Klieser W. Three-dimensional cell cultures: from molecular mechanisms to clinical applications. Am J Physiol. 1997;273:C1109–23. [PubMed]
4. Mueller-Klieser W. Tumor biology and experimental therapeutics. Crit Rev Oncol Hematol. 2000;36:123–39. [PubMed]
5. Gottfried E, Kunz-Schughart LA, Andreesen R, Kreutz M. Brave little world: spheroids as an in vitro model to study tumor-immune-cell interactions. Cell Cycle. 2006;5:691–695. [PubMed]
6. Friedrich J, Ebner R, Kunz-Schughart LA. Experimental anti-tumor therapy in 3-D: spheroids–old hat or new challenge? Int J Radiat Biol. 2007;83:849–871. [PubMed]
7. Lin RZ, Chang HY. Recent advances in three-dimensional multicellular spheroid culture for biomedical research. Biotechnology Journal. 2008;3:1172. [PubMed]
8. Hirschhaeuser F, Menne H, Dittfeld C, West J, Mueller-Klieser W, et al. Multicellular tumor spheroids: an underestimated tool is catching up again. J Biotechnol. 2010;148:3–15. [PubMed]
9. Casciari JJ, Sotirchos SV, Sutherland RM. Mathematical modelling of microenvironment and growth in EMT6/Ro multicellular tumour spheroids. Cell Prolif. 1992;25:1–22. [PubMed]
10. Venkatasubramanian R, Henson MA, Forbes NS. Incorporating energy metabolism into a growth model of multicellular tumor spheroids. J Theor Biol. 2006;242:440–53. [PubMed]
11. Sander LM, Deisboeck TS. Growth patterns of microscopic brain tumors. Phys Rev E Stat Nonlin Soft Matter Phys. 2002;66:051901. [PubMed]
12. Stein AM, Demuth T, Mobley D, Berens M, Sander LM. A mathematical model of glioblastoma tumor spheroid invasion in a three-dimensional in vitro experiment. Biophys J. 2007;92:356–65. [PubMed]
13. Jiang Y, Pjesivac-Grbovic J, Cantrell C, Freyer JP. A multiscale model for avascular tumor growth. Biophys J. 2005;89:3884–94. [PubMed]
14. Schaller G, Meyer-Hermann M. Multicellular tumor spheroid in an off-lattice Voronoi-Delaunay cell model. Phys Rev E Stat Nonlin Soft Matter Phys. 2005;71:051910. [PubMed]
15. Kim Y, Stolarska MA, Othmer HG. A Hybrid Model For Tumor Spheroid Growth In Vitro I: Theoretical Development And Early Results. Mathematical Models and Methods in Applied Sciences. 2007;17,(Suppl.):1773–1798.
16. Engelberg JA, Ropella GEP, Hunt CA. Essential operating principles for tumor spheroid growth. BMC Syst Biol. 2008;2:110. [PMC free article] [PubMed]
17. Nederman T, Norling B, Glimelius B, Carlsson J, Brunk U. Demonstration of an extracellular matrix in multicellular tumor spheroids. Cancer Res. 1984;44:3090–7. [PubMed]
18. Dorie MJ, Kallman RF, Rapacchietta DF, Van Antwerp D, Huang YR. Migration and internalization of cells and polystyrene microsphere in tumor cell spheroids. Exp Cell Res. 1982;141:201–9. [PubMed]
19. Chignola R, Foroni R, Franceschi A, Pasti M, Candiani C, et al. Heterogeneous response of individual multicellular tumour spheroids to immunotoxins and ricin toxin. Br J Cancer. 1995;72:607–14. [PMC free article] [PubMed]
20. Bortner CD, Cidlowski JA. Apoptotic volume decrease and the incredible shrinking cell. Cell Death Differ. 2002;9:1307–10. [PubMed]
21. Chignola R, Milotti E. A phenomenological approach to the simulation of metabolism and proliferation dynamics of large tumour cell populations. Phys Biol. 2005;2:8–22. [PubMed]
22. Chignola R, Del Fabbro A, Dalla Pellegrina C, Milotti E. Ab initio phenomenological simulation of the growth of large tumor cell populations. Phys Biol. 2007;4:114–33. [PubMed]
23. Chignola R, Milotti E. Numerical simulation of tumor spheroid dynamics. Physica A: Statistical Mechanics and its Applications. 2004;338:261–266.
24. Alberts B, Wilson JH, Hunt T. Molecular Biology of the Cell. New York: Garland Science, 5th edition; 2008.
25. Chignola R, Dalla Pellegrina C, Del Fabbro A, Milotti E. Thresholds, long delays and stability from generalized allosteric effect in protein networks. Physica A: Statistical and Theoretical Physics. 2006;371:463–472.
26. Milotti E, Del Fabbro A, Dalla Pellegrina C, Chignola R. Dynamics of allosteric action in multisite protein modification. Physica A: Statistical Mechanics and its Applications. 2007;379:133–150.
27. Chignola R, Schenetti A, Andrighetto G, Chiesa E, Foroni R, et al. Forecasting the growth of multicell tumour spheroids: implications for the dynamic growth of solid tumours. Cell Prolif. 2000;33:219–29. [PubMed]
28. Yu P, Mustata M, Turek JJ, French PMW, Melloch MR, et al. Holographic optical coherence imaging of tumor spheroids. Appl Phys Lett. 2003;83:575–577.
29. Milotti E, Del Fabbro A, Chignola R. Numerical integration methods for large-scale biophysical simulations. Computer Physics Communications. 2009;180:2166–2174.
30. Cornish-Bowden A. Fundamentals of enzyme kinetics. London: Butterworths; 1979.
31. O'Rourke J. Computational geometry in C. Cambridge, UK: Cambridge University Press, 2nd ed edition; 1998.
32. de Berg M, Cheon O, van Kreveld M, Overmars M. Computational Geometry: Algorithms and Applications. Berlin: Springer, 3rd ed edition; 2008.
33. Cgal, Computational Geometry Algorithms Library. Http://www.cgal.org
34. Guibas LJ, Knuth DE, Sharir M. Randomized incremental construction of Delaunay and Voronoi diagrams. Algorithmica. 1992;7:381–413.
35. Dwyer RA. Higher-dimensional Voronoi diagrams in linear expected time. Discrete and Computational Geometry. 1991;6:343–367.
36. Flekkøy EG, Coveney PV, Fabritiis GD. Foundations of dissipative particle dynamics. Phys Rev E. 2000;62:2140–2157. [PubMed]
37. Walenta S, Doetsch J, Mueller-Klieser W, Kunz-Schughart LA. Metabolic imaging in multicellular spheroids of oncogene-transfected fibroblasts. J Histochem Cytochem. 2000;48:509–22. [PubMed]
38. Mueller-Klieser W, Freyer JP, Sutherland RM. Influence of glucose and oxygen supply conditions on the oxygenation of multicellular spheroids. Br J Cancer. 1986;53:345–53. [PMC free article] [PubMed]
39. Alvarez-Pérez J, Ballesteros P, Cerdán S. Microscopic images of intraspheroidal pH by 1H magnetic resonance chemical shift imaging of pH sensitive indicators. MAGMA. 2005;18:293–301. [PubMed]
40. Acker H, Carlsson J, Holtermann G, Nederman T, Nylen T. Influence of Glucose and Buffer Capacity in the Culture Medium on Growth and pH in Spheroids of Human Thyroid Carcinoma and Human Glioma Origin. Cancer Res. 1987;47:3504–3508. [PubMed]
41. Sutherland RM, Sordat B, Bamat J, Gabbert H, Bourrat B, et al. Oxygenation and Differentiation in Multicellular Spheroids of Human Colon Carcinoma. Cancer Res. 1986;46:5320–5329. [PubMed]
42. Khaitan D, Chandna S, Arya MB, Dwarakanath BS. Establishment and characterization of multicellular spheroids from a human glioma cell line; implications for tumor therapy. J Transl Med. 2006;4:12. [PMC free article] [PubMed]
43. Deisboeck TS, Couzin ID. Collective behavior in cancer cell populations. Bioessays. 2009;31:190–7. [PubMed]
44. Wolfram Research I. Mathematica. Champaign, Illinois: Wolfram Research, Inc., version 7.0 edition; 2008.
45. Yuhas JM, Li AP, Martinez AO, Ladman AJ. A simplified method for production and growth of multicellular tumor spheroids. Cancer Res. 1977;37:3639–43. [PubMed]

Articles from PLoS ONE are provided here courtesy of Public Library of Science