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

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][2][3][4][5][6][7][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][10][11][12], or with numerical models based on rough simplifications of the biological settings such as cellular automata or other lattice-based models [13][14][15][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][10][11][12][13][14][15][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 the supporting text).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 the supporting text, which also includes important upgrades to the original model [21,22].
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 the supporting text).
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 SF1 in the supporting information 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 SF2 in the supporting text.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 mm 3 , 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 the supporting text 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 the supporting text).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 10 7 s), 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 0.1s, 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 SF3 in the supporting information).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 O(N ) time computational complexity [33][34][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 the supporting text 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 the supporting text), 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 SF4 in the supporting information).The description of the interaction forces is tuned to hold also during mitosis (see the supporting text and figure SF5).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 the supporting text.
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 O(N ) computational complexity only -and it does not grow when we include the cost of the Delaunay triangulation [35] -instead of the O(N 2 ) 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.
Several experiments [37][38][39][40][41][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 the table.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.
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.
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 reactiondiffusion scheme used in the simulation program performs correctly, and that the observed flows are not algorithmic artifacts.

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 N 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 2 N 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 SF8-21 and in the movie files S1-3 in the supporting information, 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:

Initialization
At start, the internal variables of all cells are set at approximate standard values.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 3 • 10 6 s (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 the supporting text), 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 10 7 loop variables), convergence is reasonably fast, because the initial concentration values are invariably very close to the final ones.

Cell motion
Cell motion is also regulated by differential equations and the solution uses a strategy based on a semiimplicit method (described in detail in the supporting text).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 celland 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 (2r 0 ) 2 where r 0 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.The supporting information text contains additional considerations on algorithmic complexity and on measured performance (see figure SF6 and figure SF7).
Additional processing to extract useful informations from the simulation data is performed with several standard tools, like Mathematica [44]

Supporting Information
Supporting text: provides technical details on the structure of the simulation program and includes tables ST1 to ST5.The spheroid is colored with trypan blue to mark dead cells, where the necrotic core is clearly visible.The agar contains the spheroid and helps in obtaining a better spherical shape with HeLa cells, but also stifles spheroid growth because it reduces the effective diffusion coefficients in the nourishing medium, so that it cannot be directly compared to the simulated spheroid in the second column of figure 2 (which has the same size), while it is similar to the larger spheroid in third column.projected on the plane of the section.At this stage, the necrotic core is contracting as dead cells gradually shrink, and this leads to a slow outward flow of the glucose stored in the extracellular spaces in this central region.We observe that such a behavior depends on the effective diffusion coefficient of glucose, and it disappears completely when the diffusion coefficient is high enough.This also suggests that the flow of glucose and other substances, like therapeutic drugs, is strongly dependent on the biochemistry and structure of extracellular spaces, and even small changes can lead to markedly different internal spheroid morphologies.(Right panel): individual cell velocities in the simulated spheroid.This is the same central section as in the left panel, and the velocity vectors are projected on the plane of the section.The length of each vector is proportional to the projected speed.The velocities in the viable rim show a coherent outward motion, while the velocities in the necrotic core show a rather orderly inward motion, with some vortices due to local residual cell proliferation.The region in-between is somewhat chaotic and the global structure of this plot mirrors that of the glucose flow shown on the left.The supporting information includes higher-quality versions of these figures and those of other flows.Geometry and topology of cell cluster  .Functional blocks of the C++ method that computes metabolic and extracellular variables.This part performs a loop that computes the solution of the nonlinear equations found in the implicit Euler integration step [29] (see also the supporting text).Although the number of variables can be quite large (more than 10 7 variables), convergence is fast, because the initial concentration values are invariably very close to the final ones.

Figure SF1 :
sketch of the cell phases included in the simulation program.

Figure SF2 :
sketch of the metabolic network.

Figure SF3 :
the geometry and topology of diffusion.

Figure SF4 :
pictorial representation of the interaction force between two cells.

Figure SF5 :
the geometry of mitosis.

Figure SF6 :
CPU time needed to simulate 1 hour, vs. the number of cells in the spheroid.

Figure SF18 :
lactate concentration .Figures SF19-21: velocity in the plane of the slice.Movie S1: development of the necrotic core.

Figure 2 .Figure 3 .
Figure 2. Snapshots of one simulated spheroid taken at different times.As the spheroid grows, a necrotic core develops in its central region, just as it happens in real spheroids.The size of the necrotic core and of the viable cell rim match real measurements.

Figure 8 .
Figure 8. Two views of the microstructure of a simulated spheroid, with about 500µm diameter and 296264 cells (183893 live cells + 112371 dead cells).(Left panel): flow of extracellular glucose along a central section of the tumor spheroid (yellow arrows) superposed on the plot of glucose concentration.The length of the arrows is proportional to the glucose flow intensityprojected on the plane of the section.At this stage, the necrotic core is contracting as dead cells gradually shrink, and this leads to a slow outward flow of the glucose stored in the extracellular spaces in this central region.We observe that such a behavior depends on the effective diffusion coefficient of glucose, and it disappears completely when the diffusion coefficient is high enough.This also suggests that the flow of glucose and other substances, like therapeutic drugs, is strongly dependent on the biochemistry and structure of extracellular spaces, and even small changes can lead to markedly different internal spheroid morphologies.(Right panel): individual cell velocities in the simulated spheroid.This is the same central section as in the left panel, and the velocity vectors are projected on the plane of the section.The length of each vector is proportional to the projected speed.The velocities in the viable rim show a coherent outward motion, while the velocities in the necrotic core show a rather orderly inward motion, with some vortices due to local residual cell proliferation.The region in-between is somewhat chaotic and the global structure of this plot mirrors that of the glucose flow shown on the left.The supporting information includes higher-quality versions of these figures and those of other flows.

Figure 9 .
Figure 9. Functional blocks of the simulation program.Program initialization is followed by a loop that performs biochemical and biomechanical calculations.This is followed by a check of the status of individual cells -this is where we decide whether a cell advances in the cell cycle, undergoes mitosis, or dies.Next the program computes the geometry and the topology of the cell cluster, and finally it outputs intermediate statistics and results.The loop continues until a user-defined stop condition is met.Some parts of the program can proceed in parallel (like metabolism and cell motion), and we can use multithreaded code.

Figure 10
Figure 10.Functional blocks of the C++ method that computes metabolic and extracellular variables.This part performs a loop that computes the solution of the nonlinear equations found in the implicit Euler integration step[29] (see also the supporting text).Although the number of variables can be quite large (more than 10 7 variables), convergence is fast, because the initial concentration values are invariably very close to the final ones.