# Inferring Growth Control Mechanisms in Growing Multi-cellular Spheroids of NSCLC Cells from Spatial-Temporal Image Data

Nick Jagiella, Benedikt Müller, Margareta Müller, Irene E. Vignon-Clementel, Dirk Drasdo

► To cite this version:

Nick Jagiella, Benedikt Müller, Margareta Müller, Irene E. Vignon-Clementel, Dirk Drasdo. Inferring Growth Control Mechanisms in Growing Multi-cellular Spheroids of NSCLC Cells from Spatial-Temporal Image Data. PLoS Computational Biology, 2016, 12 (2), pp.e1004412. 10.1371/journal.pcbi.1004412. hal-01244593v2

HAL Id: hal-01244593

https://hal.sorbonne-universite.fr/hal-01244593v2

Submitted on 20 Jun 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

CrossMark

click for updates

RESEARCH ARTICLE

# Inferring Growth Control Mechanisms in Growing Multi-cellular Spheroids of NSCLC Cells from Spatial-Temporal Image Data

Nick Jagiella<sup>1,2,3*</sup>, Benedikt Müller<sup>4</sup>, Margareta Müller<sup>5</sup>, Irene E. Vignon-Clementel<sup>2,6</sup>, Dirk Drasdo<sup>2,3,6*</sup>

1 Institute for Computational Biology, Helmholtz Zentrum München, Neuherberg, Germany, 2 INRIA Paris, Centre de recherche Inria de Paris, Paris, France, 3 Interdisciplinary Centre for Bioinformatics, Leipzig University, Leipzig, Germany, 4 Institute for Pathology Heidelberg (iPH), Heidelberg University Hospital, Heidelberg, Germany, 5 Faculty of Medical and Life Sciences, Furtwangen University, Furtwangen, Germany, 6 Laboratoire Jacques Louis Lions, Sorbonne Universités UPMC Univ. Paris 6, Paris, France

© These authors contributed equally to this work.

* dirk.drasdo@inria.fr

<figure><img src="image_3.png" /></figure>

OPEN ACCESS

**Citation:** Jagiella N, Müller B, Müller M, Vignon-Clementel IE, Drasdo D (2016) Inferring Growth Control Mechanisms in Growing Multi-cellular Spheroids of NSCLC Cells from Spatial-Temporal Image Data. PLoS Comput Biol 12(2): e1004412. doi:10.1371/journal.pcbi.1004412

**Editor:** Helen Byrne, University of Oxford, UNITED KINGDOM

**Received:** July 16, 2014

**Accepted:** June 24, 2015

**Published:** February 11, 2016

**Copyright:** © 2016 Jagiella et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** All relevant data are within the paper and its Supporting Information files, but due to large amount and file size of mirograph images some can only be provided as download under: http://ms.izbi.uni-leipzig.de/downloads/TumorSpheroidData.zip.

**Funding:** The authors gratefully acknowledge support by DEASE (EU FP6, No: MEST-CT-2005-021122), LungSys (BMBF, No 031602H), INVADE (INSERM, Plan Cancer, Call 2012), and NOTOX (EU-FP7). The funders had no role in study design, data

Abstract

We develop a quantitative single cell-based mathematical model for multi-cellular tumor spheroids (MCTS) of SK-MES-1 cells, a non-small cell lung cancer (NSCLC) cell line, growing under various nutrient conditions: we confront the simulations performed with this model with data on the growth kinetics and spatial labeling patterns for cell proliferation, extracellular matrix (ECM), cell distribution and cell death. We start with a simple model capturing part of the experimental observations. We then show, by performing a sensitivity analysis at each development stage of the model that its complexity needs to be stepwise increased to account for further experimental growth conditions. We thus ultimately arrive at a model that mimics the MCTS growth under multiple conditions to a great extent. Interestingly, the final model, is a minimal model capable of explaining all data simultaneously in the sense, that the number of mechanisms it contains is sufficient to explain the data and missing out any of its mechanisms did not permit fit between all data and the model within physiological parameter ranges. Nevertheless, compared to earlier models it is quite complex i.e., it includes a wide range of mechanisms discussed in biological literature. In this model, the cells lacking oxygen switch from aerobe to anaerobe glycolysis and produce lactate. Too high concentrations of lactate or too low concentrations of ATP promote cell death. Only if the extracellular matrix density overcomes a certain threshold, cells are able to enter the cell cycle. Dying cells produce a diffusive growth inhibitor. Missing out the spatial information would not permit to infer the mechanisms at work. Our findings suggest that this iterative data integration together with intermediate model sensitivity analysis at each model development stage, provide a promising strategy to infer predictive yet minimal (in the above sense) quantitative models of tumor growth, as prospectively of other tissue organization processes. Importantly, calibrating the model with two nutriment-rich growth conditions, the outcome for two nutriment-poor growth conditions could be predicted. As the final model is however quite complex, incorporating many mechanisms, space, time, and stochastic processes,
collection and analysis, decision to publish, or
preparation of the manuscript.

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

parameter identification is a challenge. This calls for more efficient strategies of imaging and
image analysis, as well as of parameter identification in stochastic agent-based simulations.

Author Summary

We here present how to parameterize a mathematical agent-based model of growing MCTS almost completely from experimental data. MCTS show a similar establishment of pathophysiological gradients and concentric arrangement of heterogeneous cell populations as found in avascular tumor nodules. We build a process chain of imaging, image processing and analysis, and mathematical modeling. In this model, each individual cell is represented by an agent populating one site of a three dimensional un-structured lattice. The spatio-temporal multi-cellular behavior, including migration, growth, division, death of each cell, is considered by a stochastic process, simulated numerically by the Gillespie algorithm. Processes on the molecular scale are described by deterministic partial differential equations for molecular concentrations, coupled to intracellular and cellular decision processes. The parameters of the multi-scale model are inferred from comparisons to the growth kinetics and from image analysis of spheroid cryosections stained for cell death, proliferation and collagen IV. Our final model assumes ATP to be the critical resource that cells try to keep constant over a wide range of oxygen and glucose medium concentrations, by switching between aerobic and anaerobic metabolism. Besides ATP, lactate is shown to be a possible explanation for the control of the necrotic core size. Direct confrontation of the model simulation results with image data on the spatial profiles of cell proliferation, ECM distribution and cell death, indicates that in addition, the effects of ECM and waste factors have to be added to explain the data. Hence the model is a tool to identify likely mechanisms at work that may subsequently be studied experimentally, proposing a model-guided experimental strategy.

Introduction

In early development, tumors grow up to 1–2mm in diameter, nourished by the nutrients and oxygen provided by the existing vasculature. Either 2D or 3D cell culture systems are utilized as biological models to study that phase, or aspects usually occurring in later phases of tumor growth and development. Current 2D cell culture approaches are only of limited use to investigate tumor progression in these stages, as they neglect crucial histo-morphological and functional features of these avascular micro-metastases or inter-capillary micro-regions of solid in vivo tumors. During the last decades, great effort has been undertaken to generate biological 3D models that describe the early phases of tumor development in a tissue context more accurately. They can thus serve as intermediate systems between traditional 2D cell culture and complex in vivo models ([3, 4]). Of these approaches, Multicellular Tumor Spheroids (MCTS) offer easy handling and fast generation, even for larger batches, and automation ([5, 6]). MCTS as a model system can be well characterized and have been shown to reproduce the spatial organization and micro-environmental factors of in vivo micro tumors, such as relevant gradients of nutrients and other molecular agents and deposition of Extracellular Matrix (ECM) (see Fig 1) ([7, 8]). Furthermore, gene expression studies revealed substantial differences in both the baseline profiles and profiles after stimulation between 2D and MCTS cultures. The latter is decisively closer to patients profile gene expression ([9–11]). Consequently, MCTS
<figure><img src="image_4.png" /><figcaption>**Fig 1. Radial organization of spheroids.** Combination of images of a cross-section of a tumor simulated with the model presented in this article depicting the spatial organization of cellular phenotypes (proliferating, dying) and the molecular agents considered by the proposed model as main resources (oxygen, glucose), growth/viability promoters (GP/VP) or growth/viability inhibitors (GI/VI). The arrows point into the direction from high to low concentrations. The image shall be compared with the corresponding scheme in ref. [8], which shows a combination of images of spheroid median sections studied with different technologies: autoradiography, TUNEL assay, bioluminescence imaging, and probing with oxygen micro-electrodes.</figcaption></figure>

doi:10.1371/journal.pcbi.1004412.g001

have now been established as experimental systems for both basic research and high through-
put screening of clinically relevant drugs (reviewed by: [12, 13]). As already explained in [14]
the organization of the different cell phenotypes within a MCTS (growing, quiescent and dead)
is supposed to be radial, and to be controlled by different factors: growth promoters (GP), via-
bility promoters (VP), growth inhibitors (GI) and viability inhibitors (VI). In the case of spher-
oids the promoters are mainly delivered from the growth medium surrounding the tumor
(with exception of ECM) while inhibitors are generally assumed to be produced by the tumor
itself. As a consequence, the local composition and interplay of those factors favor different
transitions between cell phenotypes at different distances from the tumor border.

To understand the dynamics of avascular tumor growth, several mathematical models were proposed linking the growth kinetics on the multi-cellular level (radius/volume in time) with mechanisms on cell or subcellular scales (cell growth, contact-inhibition, nutrient limitation etc.). They can be classified in two main approaches: (1) continuum models of the different components densities that evolve in time and space following PDEs (see e.g. [15–18]), (2) agent-based models that describe each cell individually and how it grows, divides, moves and dies (see e.g. [19–21]). When cells are modeled as agents, oxygen, nutrients and/or growth factors or inhibitors are often modeled by continuum models [1, 2, 22–31] or with simplified assumed profiles [32, 33]. Hybrid models on the other hand combine within the same framework the two model types for the cells, depending on the tumor zone (e.g. [26]).
Despite the large variety of models, identification of a plausible mechanistic model and its quantitative parameterization able to quantitatively explain a large set of data and predict the outcome of experiments that were not used to calibrate the model remains a difficult task. Issues are the large number of parameters and the lack of validation of the underlying mechanisms. Different models have successfully been fitted to the growth dynamics of cell populations (e.g. [2, 22] consider only the population size but not the diameter, [1, 34] consider both) relying on different mechanisms but leading to the same growth curves. For example the transition from exponential to linear radius growth phase can be due to contact-inhibition or nutrient-limitation. Based on the growth curves alone, model selection could not be made, indicating that development of a mathematical model only relying on growth curves is insufficient. In this paper we pursue a quantitative image-based approach based on bright field micrographs. This is in-line with a recent trend in large-scale simulations of brain tumors based on magnetic resonance imaging (MRI) [35, 36] following inspiring work by Swanson and coworkers [37]. Histological information has been used recently by Frieboes et. al. [38], who included histological staining measurements in a partial-differential-equation tissue model of Non-Hodgkin lymphoma growth, and by Macklin et. al. [39–41] who developed a multiscale model, mimicking cells as individual agents subject to forces, to predict ductal carcinoma growth in individual patients, including histological information.

In this paper we study mathematical model development and model parameterization by comparison with experimental data for different oxygen and glucose concentrations in the medium for the non-small cell lung cancer (NSCLC) cell line SK-MES-1. These data consist in the growth kinetics and the corresponding spatial staining patterns for nuclei, different cell states and cell environment, namely, HOECHST for cell nuclei, KI67 for proliferation, TUNEL for dead cells, collagen IV for ECM. We study in how far mechanisms that have not been directly assessed can be inferred by simultaneous matching of simulation results with experimental results on many experimental observables. Our strategy is stepwise: we first develop a model for one growth condition only, and then expand the model to capture additional growth conditions after verifying that the previous (simpler) model stage was incapable of explaining the added growth data. For this, we perform many computer simulations with the “previous” model varying each model parameter within its physiological range. Finally we arrive at a model that can almost completely be projected to the experimentally derived scheme on spheroid growth (compare Fig 1 to reference [8]). By such a stepwise strategy involving experiments, imaging, image analysis and modeling, an order mechanism during liver regeneration could be identified [42], indicating that such a strategy may be powerful in unveiling interplaying mechanisms in multi-cellular organization.

# Results

## Experiments and Image Processing

In order to study the influence of environmental conditions on the growth dynamics of tumor spheroids (Fig 2), SK-MES-1 cells were cultivated in-vitro as multi-cellular tumor spheroids under different nutriment conditions, with the hanging drop method. Then, at different points in time, spheroid size was determined with bright field microscopy and some of the spheroids were frozen, cryosectioned, stained and imaged with fluorescence microscopy (Fig 3). For a detailed description see the Materials and Methods section.

**Growth curves.** Fig 2 shows the temporal growth curves for the spheroid radius for four different culture conditions (see Table 1). The figure suggests the existence of different growth phases. The growth phase in the first two days cannot clearly be characterized but concluding from earlier work ([1, 43]) it may emerge from a transition of an exponential to a sub-
<figure><img src="image_5.png" /><figcaption>**Fig 2. Growth curves of MCTS cultivated under different nutrient conditions.** The arrows and boxes indicate the time points where histological images were taken.</figcaption></figure>

doi:10.1371/journal.pcbi.1004412.g002

exponential growth phase. Between days 2–21, growth is over all approximately linear, followed by saturation or even shrinkage of the spheroid. The crossover points between the growth phases depend on the oxygen and glucose medium concentrations. Hence there must be an inflection point between the initial convex and the final concave (saturation) phase where the 2nd derivative of the radius with respect to time changes from positive to negative values. For [G] = 5mM, 25mM, [O] = 0.28mM, the growth behavior is very similar. Lowering glucose to [G] = 1mM or oxygen to [O] = 0.07mM leads to significantly slower growth. At about 24 days the radius is saturating or even shrinking.

**Radial profiles of histological quantification.** In order to obtain spatial information on the distribution of cells, proliferating cells, dead cells and extracellular matrix (ECM), several multi-cellular spheroids were stained and imaged. Fig 3 shows typical images stained with HOECHEST, Ki67, TUNEL and Collagen IV for cell nuclei, proliferating cells, dying cells and extra-cellular matrix (ECM), respectively for [G] = 5mM and [G] = 25mM. The images were segmented and cells classified into being proliferating (if the nucleus was Ki67—positive), dying (if TUNEL—positive) or quiescent (if neither Ki67 nor TUNEL positive) (Fig 4, details see in Materials and Methods). From the segmented images the spatial profiles for proliferating, quiescent and dead cells were investigated, from the border to the interior (for details see section Quantitative Image Analysis). Starting the investigation of spatial profiles from the border (as opposed to from the center of mass) ensured that the profiles at the border were correctly captured even if the MCTS had an ellipsoidal instead of a spherical shape, which was often the case (see for example, Fig 4A and 4E). The disadvantage of this choice was that the number of events (Ki67 positive nuclei, ECM density etc.) inside the tumor was usually very rare so that the data points with distance of more than 200 micro-meters from the border have to be considered with great caution (far inside the tumor the number of events is very small and the fraction depends on the bin-size, see Fig S13 in S1 Document). The fraction of proliferating cells decreases the farther the proliferating cells are inside the spheroid (Fig 3). The proliferation profiles do (almost)
<figure><img src="image_6.png" /><figcaption>**Fig 3.** Quantification of proliferating (top) and dying (bottom) cell nuclei as well as of extra-cellular matrix density (center). The images (top) depict cross-sections taken at different time points (T3 = 17d, T4 = 24d) of spheroids grown under different nutrient conditions ([G], [O₂]). The colors indicate cell nuclei (blue), extra-cellular matrix (EMC) (green) and proliferation (red, left) or cell death (red, right). The curves (bottom) represent the radial profiles of the fraction of proliferating and dying cells, and the ECM density (intensity of Collagen IV staining), inferred from averages over images from several spheroids growing under the same conditions. Bars indicate the standard deviation. With ECM density we denote the fluorescence value of Collagen IV, which varied in the interval [0, 1].</figcaption></figure>

doi:10.1371/journal.pcbi.1004412.g003

not change over time and are very similar for [G] = 5mM and [G] = 25mM in the growth medium. Cell death increases from the spheroid border to its interior. There is a marked difference in the cell death profile between both nutrient conditions in that for 5mM the necrotic core emerges earlier (already at day 17, while for 25mM not until day 24). Interestingly after 24 days,
<table><caption>Table 1. Nutrient conditions.</caption><thead><tr><th>Condition</th><th>[G]</th><th>[O<sub>2</sub>]</th><th>Note</th></tr></thead><tbody><tr><td>I</td><td>1mM</td><td>0.28mM</td><td>hyponourished</td></tr><tr><td>II</td><td>5mM</td><td>0.28mM</td><td>physiological</td></tr><tr><td>III</td><td>25mM</td><td>0.28mM</td><td>hypernourished</td></tr><tr><td>IV</td><td>25mM</td><td>0.07mM</td><td>hypoxic</td></tr></tbody></table>

doi:10.1371/journal.pcbi.1004412.t001

cell death seems to occur at smaller penetration depth for [G] = 25mM than for [G] = 5mM. The TUNEL assay does not permit distinguishing between necrotic and apoptotic cells ([44]). For this reason we also did some pilot experiments selectively staining samples with Caspase 3. In these samples we could see activation of Caspase 3 in the peripheral parts of what we would identify as the “necrotic core”, although this may also not yield a definite distinction between apoptosis and necrosis, as Caspase 3 activation has been shown to be involved in both ([45, 46]). Thus we could not discriminate between different forms of cell deaths in the experiments within this paper. Indeed, as this cell culture system does not include macrophages, removal of apoptotic cell bodies is very slow so that distinguishing between necrotic and apoptotic cells seemed dispensable. This assumption was self-consistently confirmed as very small apoptosis and lysis rates had to be assumed in the model to explain the experimental findings so that apoptosis did not have to be taken explicitly into account.

The amount of extra-cellular matrix (ECM) increases with the penetration depth until 100–200µm. For larger penetration depths it seems to decrease again, but the data is too noisy to infer a clear course. Counterintuitively, cells at the spheroid surface have a significantly lower proliferative activity than those a few layers inside.

## Mathematical Modeling: Main Components Needed

Our objective is to explain the experimentally observed growth pattern for different glucose and oxygen medium concentrations within one mathematical model. For this purpose we first searched for a minimal model (in the sense specified in the abstract) explaining the experimental tumor growth observations for medium concentrations of [G] = 25mM and [O] = 0.28mM, and then stepwise extended this model to capture the other growth conditions (for illustration of the stepwise model development strategy, see Fig S3 in S1 Document). We based our choice of possible control mechanisms upon prior knowledge guided by published information and own experiments. We have chosen the condition of maximum glucose and oxygen medium concentration as the similarity of the growth kinetics for [G] = 25mM and [O] = 28mM versus [G] = 5mM and [O] = 28mM suggests that for the former condition neither glucose nor oxygen may be limiting (see also Discussion below). This line of argument is supported by findings of Freyer and Sutherland for another cell type at almost the same oxygen and glucose medium conditions (compare with [1]). We fit at each model development step all parameters again. So the fits shown in this article were the best we could obtain for the respective model. However, due to the large search space and duration of simulation of at least one day (reference computer: Intel(R) Xeon(R) processor X5680 3.33GHz 12M cache 6-core and 144 GB DDR3-RAM 1333 MHz) it cannot completely be excluded that further parameter searches could give additional slight improvements. In order to promote readability, we enumerate the model at each development level.

Usually we performed for each parameter set only a single simulation i.e., a single realization of the stochastic growth process. This can be justified by observing that the growth process, that starts with about 10000 cells as in the experiment, is self-averaging such that the variations
<figure><img src="image_7.png" /><figcaption>**Fig 4. Image smoothing and segmentation.** (A) The micrograph shows a cryosection of a spheroid stained with HOECHST (red), Ki67 (blue) and Collagen IV (green). The original image was smoothed with a median filter. For a zoomed-in section of the image (B), the original (C) and smoothed version (D) of the image are visualized as landscapes of the inverted blue color channel intensities, $1 - I^{blue}$. (E) The cell nuclei were segmented from the blue color channel and differentiated between Ki67 positive (red) and negative (blue) nuclei by means of the red channel. The spheroid lumen (green) is approximated by inflated nuclei.</figcaption></figure>

doi:10.1371/journal.pcbi.1004412.g004
<figure><img src="image_8.png" /></figure>

for different realizations of the growth process for the same parameters are negligible (Fig S11 in S1 Document).

Our basic model considers each cell individually within an agent-based model i.e. each individual cell is represented by an agent. Molecules, which finally are glucose, oxygen and lactate, as well as extra-cellular matrix and “waste” material released by dying cells, are represented by their local concentrations. We use the term “molecules” in what follows to generically describe these environmental factors that affect the cells. The model is three-dimensional. In the following sections we will introduce briefly the main model components. A detailed description of the model as well as the biological processes mimicked can be found in the material and method section and in the supporting information (S1 Document).

**Cells.** The cells were modeled as individual objects populating the sites of an unstructured lattice (see supplement, [25, 47]) with no more than one cell on a lattice site. The cells could be in either of three states: proliferating, quiescent or dying. Depending on their state and their local environment they were able to progress in or reenter (from G0) into the cell cycle, grow, migrate, die or undergo lysis (Fig 5). Each of these processes required a set of necessary conditions, and was executed with a certain rate. Both, processes and rates depended on the cellular neighborhood and the local molecular concentrations (see below). The temporal dynamics of the cells was described by a “master equation” mimicking the time evolution of the probability

<figure><img src="image_9.png" /><figcaption>**Fig 5.** Scheme of cell states and the biological processes a cell can undergo: cell cycle (including cell growth and division), migration, cell death (apoptosis and necrosis) and lysis. Left panel: A cycling cell grows by occupying two neighboring lattice sites after a fraction of transitions in the cell cycle. Growing cells can push a certain number of cells aside. Cells can migrate by hopping from one lattice site to its neighbor lattice site. A cell dies with rate $k_{\text{neck}}$ and is consequently lysed with rate $k_{\text{lys}}$. (The schemes in the left panel are 2D for clarity; the simulations however were all in 3D). Right: All processes are modeled as Poisson processes. A cell in the cell cycle undergoes $m_d$ transitions until splitting into two daughter cells. Composed of $m_d$ sub-processes the cell cycle time ends up following an Erlang distribution. Increasing $m_d$ will lead to sharper cell cycle distributions. A cell transitions from one cell cycle state (CCS) to the next with rate $k_{\text{div},m}$. If a cell is in state $m_g$ it will grow. If it is in state $m_d$ it will divide into 2 daughter cells during the next transition. Furthermore, the 2 daughter cells will either enter the first CCS proportional to probability $p_{\text{div}}$ or become quiescent (G0) proportional to probability $1 - p_{\text{div}}$. A cell in quiescence can reenter the cell cycle with rate $k_{\text{re}}$ and probability $p_{\text{re}}$. Table 2 indicates how all transition rates and probabilities are calculated for models 1–4.</figcaption></figure>

doi:10.1371/journal.pcbi.1004412.g005
of a certain multicellular configuration. This equation was numerically simulated with the Gillespie algorithm (see S1 Document).

**Molecules.** For each molecule type a continuum representation of its concentration was chosen. Its spatial-temporal dynamics was modeled by a partial differential equation (PDE). For each molecule the equation could be casted into the form of a reaction diffusion equation

$$ \partial_t u = \nabla \cdot (D_u(\sigma)\nabla u) + r_u(\sigma, u), \quad (1) $$

where $t$ is the time, $u \in \{G, O, L, ECM, W\}$ is the concentration of either glucose, oxygen, lactate, extra-cellular matrix or waste, $D_u$ is the respective diffusion coefficient, and $r_u$ a rate function representing the reactions involving $u(= u(\mathbf{r}, t))$. $\mathbf{r}$ denotes the position in space. The diffusion coefficients as well as the reaction functions depended on the local cell density $\sigma(\mathbf{r}, t) = \sum_k \delta(\mathbf{r} - \mathbf{r}_k)$, $\mathbf{r}_k(t)$ being the position of cell $k$ at time $t$. The reactions were specific to each component (molecules or ECM): cells consumed glucose and oxygen with $-r_G$ and $-r_O$, lactate was produced by under-oxygenated cells with $r_L$, ECM was produced by cells and disintegrated with a certain rate, and waste was released by dying cells.

## Mathematical Modeling: Coupling of Scales and Comparison with Data

A number of parameters related uniquely to cells or molecules were either taken or estimated from literature (e.g. molecular diffusion coefficients, consumption rates, cell cycle time distribution). Others could be inferred from the data presented in this article either directly (e.g. initial conditions, cell size) or by sensitivity analysis (e.g. cellular division, cell death and lysis rates) (see Table S1 in S1 Document). The identification of the main mechanisms coupling the cellular and molecular kinetics and their parameterization (see Fig 5) was subject to the model comparison with data explained below.

The labeling patterns in Figs 3 and 4 confirm a border distance-dependent “zonation” as discussed in the introduction (compare also [8] and Fig 1). In the following we will study the influence of cell-cell-contacts, extra-cellular matrix and metabolic compounds (nutrients/metabolites) to infer the corresponding model mechanisms. For the latter (metabolism) we will compare four different hypotheses (model 1–4). A summary of the equations used to calculate the transition rates and probabilities depicted in Fig 5 for all models studied in the following sections can be found in Table 2.

<table><caption>**Table 2. Transition rates and probabilities of models 1–4.** All models have $p_{re} k_{re} = 0$, $k_{apt} = \text{const.}$ and $k_{lys} = \text{const.}$ in common. The transition rates $k_{div,m}$ between individual cell cycle steps is given by $k_{div,m} = m_d k_{div}$. All parameter values can be found in Table S1 in S1 Document.</caption><thead><tr><th>Model</th><th><i>p</i><sub>div</sub> =</th><th><i>k</i><sub>div</sub> =</th><th><i>k</i><sub> nec </sub> =</th></tr></thead><tbody><tr><td>1</td><td>e<sup>-L/ΔL</sup></td><td>k<sup>max</sup><sub> div </sub></td><td>k<sup>max</sup><sub> nec </sub></td></tr><tr><td></td><td>·H([ECM] - [ECM]<sup>min</sup>)</td><td>·H([G][O] - p<sup>oxyglc</sup>)</td><td>·H(p<sup>oxyglc</sup> - [G][O])</td></tr><tr><td>2</td><td>e<sup>-L/ΔL</sup></td><td>k<sup>max</sup><sub> div </sub></td><td>k<sup>max</sup><sub> nec </sub></td></tr><tr><td></td><td>·H([ECM] - [ECM]<sup>min</sup>)</td><td>·H(p<sub>ATP</sub> - p<sup>min</sup><sub> ATP </sub>)</td><td>·H(p<sup>min</sup><sub> ATP </sub> - p<sub>ATP </sub>)</td></tr><tr><td>3</td><td>e<sup>-L/ΔL</sup></td><td>k<sup>max</sup><sub> div </sub></td><td>k<sup>max</sup><sub> nec </sub></td></tr><tr><td></td><td>·H([ECM] - [ECM]<sup>min</sup>)</td><td>·H(p<sub>ATP</sub> - p<sup>min</sup><sub> ATP </sub>)</td><td>·H(p<sup>min</sup><sub> ATP </sub> - p<sub>ATP </sub>)</td></tr><tr><td>4</td><td>e<sup>-L/ΔL</sup></td><td>k<sup>max</sup><sub> div </sub></td><td>k<sup>max</sup><sub> nec </sub></td></tr><tr><td></td><td>·H([ECM] - [ECM]<sup>min</sup>)</td><td>·H(p<sub>ATP</sub> - p<sup>min</sup><sub> ATP </sub>)</td><td>·H(p<sup>min</sup><sub> ATP </sub> - p<sub>ATP </sub>)</td></tr><tr><td></td><td>·H(n<sup>max</sup><sub> exp </sub> - n<sup>exp</sup>)</td><td>·H(1 - 0.5([W] - [W]<sup>max</sup>))</td><td>·H(1 - 0.5([O]<sup>min</sup> - [O]))</td></tr></tbody></table>

doi:10.1371/journal.pcbi.1004412.t002
**Mechanical growth-inhibition: Probabilistic cell decision.** The growth curves and proliferation profiles for [G] = 25mM and [G] = 5mM look quite similar while the necrotic core sizes at 17 days look very different for [G] = 25mM and [G] = 5mM (Figs 2 and 3). An equivalent observation has been made for data on EMT6/Ro cells [14] and could be shown by model simulations to strongly indicate that proliferation is not controlled by non-glucose-related factors but by a mechanical form of contact inhibition [1]. Otherwise, the proliferating rim and the size of the necrotic zone should both be affected. Model simulations with our lattice model confirmed this finding.

In previous articles a biomechanical form of contact-inhibition was assumed to be a deterministic decision depending on the local pressure or deformation in growing cell populations ([1, 2, 48–50]). Nevertheless, as can be seen from model studies of growing multi-cellular spheroids in an off-lattice model in [24] and [34], a deterministic criterion based on a deformation threshold is able to yield a gradual change of the fraction of dividing cells as a consequence of local stochastic events such as cell micro-motility or randomness in the cell cycle duration, which both introduce stochastic fluctuations in deformation and mechanical stress. As demonstrated in refs. ([25, 43]) an extended proliferating rim can be qualitatively captured within a lattice model (as the one chosen in this article) by introducing a maximum distance, ΔL, over which a cell is able to push other cells away in order to grow. This probability that a cell enters the cell cycle can then be defined as follows

$$p_{\text{div}} := H(L - \Delta L), \qquad (2)$$

where H is the Heaviside function, H(x) = if x < 0, H(x) = 1 if x ≥ 0, and L is the distance to the closest free lattice site.

However, as illustrated in Fig 6 (left panel) the introduction of such a sharp cut-off leads to a relatively sharp transition between outer proliferating and inner quiescent rim while the experimental profiles display a wide transition zone with gradual change of the frequency of

<figure><img src="image_10.png" /><figcaption>Fig 6. Experimental and model-predicted fraction of Ki67 positive cells (left) and Collagen IV (ECM) density (right) for [G] = 25mM, [O] = 0.28mM at day 17 (red) and day 24 (green) vs. distance from the spheroid surface. A deterministic cell cycle progression depending on the distance of a cell from the surface, generates sharp transitions that are experimentally not observed (left, dotted line), while a probabilistic transition generates a smooth monotonically decreasing function of Ki67 positive cells in the simulations (left, dashed). If in addition cells can progress in the cell cycle only if the local concentration of ECM overcomes a threshold value [ECM]min, the experimentally found Ki67 profile is recovered (left, full line). Right full line: ECM concentration in the computer simulation.</figcaption></figure>

doi:10.1371/journal.pcbi.1004412.g006
cells in the cell cycle. The smooth experimentally observed transition could be obtained within a lattice model by choosing a probabilistic approach instead, with cells entering the cell cycle with a probability exponentially decreasing with $\Delta L$ (Fig 6, left panel):

$$p_{\text{div}} := e^{-L/\Delta L} \qquad (3)$$

Hence the decision of a cell entering the cell cycle is stochastic, meaning that neighbor cells with the same distance to the border may behave differently. This generates local inhomogeneous cell behavior. The cell cycle progression probability (3) generates a profile with the highest Ki67 positive cell fraction at the outermost border. This tendency is a robust phenomenon and is also found in off-lattice models of growing cell population where the cell cycle progression is controlled by a deformation or a pressure threshold ([48, 51]). Hence we conclude that biomechanical control of cell cycle progression by a deformation or pressure threshold alone cannot explain a drop of proliferation at the margin of the multicellular spheroid.

**Extra-cellular matrix.** The experimental data exhibit that border cells have a lower proliferating activity than interior cells (Fig 6, left panel) while at the same time collagen IV staining display a lower ECM concentration close to the spheroid border (Fig 6; right panel). This suggests that the known control function of ECM in the cell cycle progression [50] may be responsible for the lower fraction of proliferating cells at the spheroid border. This assumption is supported by the observation that ECM binding and connected integrin signaling is an important part of proliferative signaling ([52, 53]).

Accordingly we assume that a minimal density of extra-cellular matrix, is required for cells to enter the cell cycle by replacing Eq 3 by

$$p_{\text{div}} := e^{-L/\Delta L} \cdot H([\text{ECM}] - [\text{ECM}]^{\text{min}}), \qquad (4)$$

where [ECM] reflects the local concentration of ECM, and [ECM]^{min} is a threshold value (the values do not have a dimension as we directly work with the fluorescence intensity defined in the interval [0, 1]). Extra-cellular matrix is secreted by each cell, and decreases with a certain rate. Choosing [ECM]^{min} = 0.003, the extracellular matrix concentration behaves as experimentally observed (Fig 6), and the simulated proliferation profile shows good agreement with the experimental profile.

Eq 4 determines the probability of cell cycle entrance of each daughter cell directly after division. For most of the paper we assumed that this decision is irreversible ($p_{re} = 0$ in Fig 5) but we also considered the consequence of exit from G0 ($p_{re} > 0$ in Fig 5). For the latter we tested different alternative rules (see Fig S4 and S5 in S1 Document) that may be subsumed by the equation for the exit rate, $r(L) = k_{re} p_{re}(L)$. This equation assumes that the cell probes with rate $k_{re}$ exit of G0 occurring with probability $p_{re}(L)$. The probability that the entrance occurs within time interval $[t, t + \Delta t)$ at depth $L(t)$ is then

$$p(\text{re - entrance in } [t, t + \Delta t)) = k_{re} p_{re} e^{-\int_{t}^{t+\Delta t} p_{re} k_{re} dt'} \Delta t.$$ In case $p_{re}$ has the functional form as $p_{div}$ in Eq 4, the rate $k_{re}$ must be very small in order to achieve agreement between data and model (Fig S4 in S1 Document). We found $k_{re} = 0.001/2.4h$ (Fig S4 in S1 Document; given by Eq 4). I.e., only in $100/p_{re}$ days on the average a specific cell exits G0.

We furthermore tested the hypothesis that cells become quiescent with probability $1 - p_{div}$ (as in Eq 4) and can only exit quiescence if their neighbor site becomes free (Fig S5 in S1 Document) which may occur if a cell is able to sense free space. In this case we can also obtain a good fit to the growth data if the threshold ECM concentration at which re-entrance can occur, is elevated. We note, that very few cell cycle reentrances from G0 could be observed during simulations, and that the simulation results were not sensitive at all to changes in the
reentrance rates (see S1 Document). We concluded that the presence of cell—cycle re-entrance by exit of the quiescent state is not a critical parameter for MCS growth of SK-MES-1 cells. However, in case of drug therapy, where a drug may kill the outermost cells, cell-cycle re-entrance is expected to modify growth and survival of the treated cell population, and would then become a critical parameter.

In conclusion the cell-cycle re-entrance must be a very rare (negligibly rare) event in order to explain the data which is why we throughout the following neglect cell cycle re-entrance (for more details, see S1 Document).

**Model 1: Glucose- and oxygen-limited growth and survival.** We assumed deprivation from glucose and oxygen to be responsible for cell death. Schaller and Meyer-Hermann [2] proposed that cell survival is possible only if the product of glucose and oxygen overcomes a threshold $p^{oxgluc}$ while below this threshold cells undergo cell death:

$$k_{div} = k_{div}^{max} H([G][O] - p^{oxgluc}), \qquad (5)$$

$$k_{ nec } = k_{ nec }^{ max } H(p^{ oxgluc } - [G][O]), \qquad (6)$$

where $H$ is the Heaviside function, $k_{div}^{max}$ and $k_{ nec }^{max}$ are the maximal and $k_{div}$ and $k_{div}$ are the actual division and cell death rates, respectively. As no experimental measurements of glucose and oxygen consumption rates for SK-MES-1 cells were accessible, we assumed the rates to obey the same functional relationship as the oxygen and glucose consumption rates measured for EMT6/Ro cells (see Materials and Methods), permitting only the parameters to be different for SK-MES-1 cells than for EMT6/Ro cells.

Fig 7 shows that for this assumption (dashed lines) and nutrient medium condition III ([G] = 25mM, [O] = 0.28mM) the simulation results are in good agreement with the experimental growth curves and the radial profiles for ECM, proliferation and cell death, but completely fail for condition II ([G] = 5mM, [O] = 0.28mM). The size of the necrotic core and the growth-lim-itation are over-estimated.

**Model 2: ATP-limited growth and survival.** In a next step we considered cell cycle pro-gression and death to depend on the yield of adenosine triphosphate (ATP), the main energy currency of a cell [54] a cell can produce from the consumed glucose and oxygen, instead of the local product of oxygen and glucose concentration:

$$k_{div} = k_{div}^{max} H(p_{ATP} - p_{ATP}^{min}), \qquad (7)$$

$$k_{ nec } = k_{ nec }^{ max } H(p_{ATP}^{min} - p_{ATP}), \qquad (8)$$

where $p_{ATP}$ is the production rate of ATP (see Eq 28), $p_{ATP}^{min} = 900mM/h$ a threshold value demarcating the border between proliferating and dead states. We note that cells which are suf-ficiently supplied with both, oxygen and glucose, produce ATP at a rate of about 1100mM/h (see Fig 8 which seems to optimally fulfill the cells energy requirement. Interestingly, cells that moderately lack either glucose or oxygen will increase their ATP production up to 1700mM/h.

Fig 7 shows that ATP-limited growth (solid lines) is in good agreement with the tumor growth kinetics for both conditions, II and III.

The view of an important growth regulating function of ATP is supported by the observa-tions that ATP has been shown to control the cell cycle in developing retinas of chick and mouse [55], and that the most successful drugs targeting the cell cycle inhibit ATP binding to Cyclin Dependent Kinases—Enzymes (CDKs) essential for cell cycle progression ([56, 57]). Moreover, ATP has been found to upregulate proliferation of human cardiac fibroblasts ([58]).
<figure><img src="image_11.png" /><figcaption>**Fig 7.** Growth kinetics (upper panel) and spatial profiles of TUNEL, Ki67, Collagen IV for [G] = 5mM, [O] = 0.28mM (condition II, upper right picture and middle panel) and [G] = 25mM, [O] = 0.28mM (condition III, upper left picture and low panel.). The dotted curves show simulations where cell cycle progression requires the product of local glucose and oxygen concentration to exceed a critical threshold, while below this threshold cells become quiescent and die. If the product is replaced by the local concentration of ATP, the agreement improves significantly (solid lines). For dashed curves the diffusion coefficient in medium was assumed to be $D_{med} = 30D_{tum}$. The experimental growth curves are shown as boxplots (box: mean, lower & upper quartiles, horizontal dashes: minimum & maximum), and the radial profiles as composition of mean (bold line) and standard deviation (thin line).</figcaption></figure>

doi:10.1371/journal.pcbi.1004412.g007

Nevertheless, in the case of 5mM of glucose (condition II) the model predicts a much larger necrotic core than for 25mM of glucose (condition III) while the experimental curves suggest a rather comparable size. Moreover, the fraction of dying cells for condition II is much smaller than experimentally found. The size of the necrotic core changes with the diffusion coefficients of glucose and oxygen as illustrated for another simulation example shown in Fig 7 (dashed lines). Here we assumed the molecular diffusion in the medium to be 30-fold larger than in the spheroids. In this scenario the model shows a better agreement with 5mM of glucose (condition II), but rather underestimates the size of the necrotic core for 25mM of glucose (condition III). Although the size of the necrotic core predicted by model 2 changes with the medium concentration of glucose and apparently the diffusion coefficients, the experimental curves suggest that it is approximately the same for 5mM and 25mM of glucose (compare also day 24 in Fig 3).

**Model 3: Lactate-induced cell death.** Moreover, the experimental profiles for the dying cell fraction (see Fig 3) reveal several differences between both conditions II and III. Firstly, for 5mM glucose the necrotic core starts at a depth of about 150 microns at 17 days as well as at 24
<figure><img src="image_12.png" /><figcaption>**Fig 8. Nutrient-dependent metabolism.** Production rates of ATP (top left) and lactate (top right). The threshold value of the minimal ATP requirement is depicted by the black line. Percentage of anaerobically metabolized glucose (bottom left) and glycolytic ATP (bottom right).</figcaption></figure>

doi:10.1371/journal.pcbi.1004412.g008

days, while for [G] = 25mM it emerges slowly: it is almost not present at day 17 and becomes big at day 24. Secondly, the transition between necrotic and viable parts is much sharper for [G] = 5mM than for [G] = 25mM. And third, after 24 days, counter intuitively, there are even more dying cells found in the exterior parts (< 150µm) for 25mM of glucose compared to 5mM even though the latter is supposed to be more nutrient-limited.

This suggests that the cause of cell death may be different for the two conditions. Hence in a third step we studied the consequence of the additional assumption that cell death may also be triggered if the local lactate concentration overcomes a critical threshold (Fig 9). Cells being sufficiently supplied by glucose, but lacking oxygen, produce lactate. We here test the assumption that lactate might have a death-promoting effect by setting the death rate to

$$k_{\text{nec}} = k_{\text{nec}}^{\text{max}} \cdot H(p_{\text{ATP}}^{\text{min}} - p_{\text{ATP}}) \cdot \frac{[L]^{n}}{([L]^{\text{max}})^{n} - [L]^{n}}, \quad (9)$$

where $[L]^{\text{max}} = 20mM$ is the lactate concentration where the cell death rate is half its value, and $n$ is the Hill coefficient which determines the sharpness of the switch from no to full cell death response. I.e. $n = 1$ leads to a smooth transition and $n \to \infty$ to a Heaviside function. We have chosen the Hill coefficient such that the best fit could be obtained (Fig S8 in S1 Document). A linear dependency or a Heaviside function led to significantly worse agreement to the data than a Hill function with $n = 2$. As lactate is a side-product of anaerobe metabolism Eq 9 will mainly affect conditions with high glucose, but low oxygen supply (see Fig 8).
<figure><img src="image_13.png" /><figcaption>**Fig 9. Adding lactate-induced cell death.** If above a certain concentration of lactate produced during glucogenesis, cell death occurs, population kinetics and spatial temporal profiles for conditions III ([G] = 25mM, [O] = 0.28mM), and II ([G] = 5mM, [O] = 0.28mM), as well as the growth phases of condition I ([G] = 1mM, [O] = 0.28mM) and condition IV ([G] = 25mM, [O] = 0.07mM) were correctly captured. The diffusion coefficient in the medium was assumed to be $D_{med} = 30 D_{tum}$. The experimental growth curves are shown as boxplots (box: mean, lower & upper quartils, horizontal dashes: minimum & maximum), and the radial profiles as composition of mean (line) and standard deviation (error bars).</figcaption></figure>

doi:10.1371/journal.pcbi.1004412.g009

Fig 9 shows how the combination of ATP-limited and lactate-induced cell death impacts the growth curves for all four conditions and the radial profiles of condition II ([G] = 5mM, [O] = 0.28mM) and III ([G] = 25mM, [O] = 0.28mM). Now both, shapes and temporal evolutions of the radial profiles of the proliferating and dying cell fractions are in much better agreement with the data than before.

Moreover, and interestingly, the model, although only calibrated with [G] = 5mM, [G] = 25mM and [O] = 0.28mM, quantitatively predicts the growth phase for very poor glucose medium concentration ([G] = 1mM, [O] = 0.28mM) and very poor oxygen medium
concentration ([G] = 25mM, [O] = 0.07mM) without any additional parameter adjustment! This indicates that the functional dependencies on the key factors determining the growth phase of SK-MES-1 spheroids are captured.

It is interesting to notice that growth adverse effects of lactate have been observed in mouse hybridoma ([59]). High concentrations of lactate (55mM) were found to inhibit cell growth, with accompanying increase of death rates (even though the authors attribute the latter to higher osmolality). However, the overall dependency of the death rate on lactate concentration is very similar to those predicted by our model (Fig S7 in S1 Document). A similar finding has been made in CHO (Chinese Hamster Cells; [60]), immortalized cells.

**Model 4: Waste and under-oxygenation mediating quiescence.** Finally, although the initial growth phases observed for all four conditions can be well explained by ATP and lactate dependent cell dynamics, they fail to explain saturation, which occurs for all conditions after about 20 days. As for small times for all growth conditions for which [O] = 0.28mM growth is faster than linear, saturating at about 20–30 days (i.e., turning into a phase in which growth is slower than linear), there must exist a point of inflection, above which the expansion speed starts to decrease. However, as long as the nutrient concentration in the medium is maintained, the outer cells are fed by nutrients and a decrease of the expansion speed cannot be nutrient related. Hence the reduction in the expansion speed must emerge from growth inhibition originating from processes inside the tumor. Lactate can be excluded as a candidate for such a growth inhibitor as it only appears in the cases of oxygen-deprivation (III & IV). Moreover, mechanical stress can be excluded to be responsible for saturation as the environment of the multi-cellular spheroid is liquid, not a medium that responds with higher resistance upon compression ([48, 61, 62]).

We tested the hypothesis that the cytoplasmic waste material, released by lysing dying cells, may act as growth inhibitor, by the assumption that cells being exposed to waste for longer than a certain time span become quiescent:

$$n_{exp} = \int_{-\infty}^{t} [1 - H([W]^{max} - [W](\mathbf{r}, \tau)) \cdot H([O](\mathbf{r}, \tau) - [O]^{min})] \cdot k_{div} d\tau \geq n_{exp}^{max}, \quad (10)$$

where H is the Heaviside function, $n_{exp}$ is the number of cell cycles a cell was either exposed to waste material or deprived of oxygen, or both, $n_{exp}^{max}$ is the number of cell cycles a cell can maximally stay under this condition without being growth inhibited, $[W](\mathbf{r}, \tau)$ and $[O](\mathbf{r}, \tau)$ the local concentrations of waste and oxygen molecules, $[W]^{max}$ and $[O]^{min}$ the respective threshold concentrations. I.e., a cell that has accumulated too much waste becomes quiescent, and does not (re-)enter the cell cycle anymore. Notice, that this dynamics can also be expressed by a stochastic dynamics and hence incorporated into the master equation by introducing an internal state variable (counter) for each cell that increases with a certain constant rate when the external local concentration of waste overcomes the threshold concentration, and assuming cell cycle (re-)entrance to be inhibited when this state variable overcomes a critical value. If the number of states and the rate of increasing the counter by one unit are both large, the emerging waiting time until the threshold is reached is Erlang-distributed with very small standard deviation (i.e. comes very close to the deterministic limit represented by Eq 10). Instead of counting the number of cell cycles the cell is exposed to a critical concentration of waste or to insufficient oxygen (or both) as in Eq 10, we tested alternatively the condition, that the exposure is integrated over absolute time (instead of multiple of the cell cycle times). In this case, only the fit to condition IV changes by that the deviation between experiment and simulation increases by about 5%.
In summary, the transitions between proliferating and quiescence for model 4 are determined by Eqs 4 and 10.

Moreover, it can be observed that the radius expansion slows down under all four conditions after about 200...300 hours. Thus, we further assume that cells under waste exposure or oxygen deprivation grow with a reduced speed of

$$k_{\text{div}} = k_{\text{div}}^{\text{max}} \cdot H(p_{\text{ATP}} - p_{\text{ATP}}^{\text{min}}) \cdot (1 - 0.5H([W] - [W]^{\text{max}})) \cdot (1 - 0.5H([O]^{\text{min}} - [O])) \quad (11)$$

where $[O]^{\text{min}}$ and $[W]^{\text{max}}$ are the critical concentrations of oxygen and waste, respectively. The growth rate reduction for low oxygen supply was previously reported by Stolper et al. [63], while the growth-inhibiting effect of waste originating from necrotic debris was already proposed by Greenspan [64].

Fig 10 shows the simulation results with the extended model, including the lactate and waste dependency for two different thresholds of the waste exposure time. Condition II and III show a similar initial growth phase, followed by two linear regimes—the first dominated by ATP-limitation and the second by waste-intoxication—and finally terminating in saturation.

On the other hand, conditions I and IV saturate almost immediately after a very short initial growth phase driven by waste exposure (condition I) or oxygen deprivation (condition IV). Note that saturation is only intermittent if cell quiescence is delayed by a sufficient number of cell cycles under exposure. While for a larger delay of 8 cell cycles (Fig 10) the spheroids permanently stop to grow, for a 4 cell cycle delay (not shown) the spheroid expansion would only temporarily slow down and then recover its initial pace. This is related to the larger amount of waste, which can pile up and spread during longer delay times and thus affect all cells and not only the ones closest to the necrotic core (for simulation snapshots of all conditions see Fig S14 in S1 Document).

## Discussion

In this paper we inferred a mathematical model of tumor spheroid growth for the non-small cell lung cancer cell line SK-MES-1 from image data of growing tumor spheroids. Cell nuclei, proliferating cells, extra-cellular matrix and dying cells (by either necrosis or apoptosis) were labeled at different points in time and under different oxygen and glucose medium concentrations. The model was built by an iterative procedure, which we propose as a general template for modeling tissue organization processes. We started by developing a minimal model for one growth condition only, then stepwise extending this model by further mechanisms whenever the previous simpler model turned out to be insufficient to reproduce the experimental observations for an additional growth condition. Before adding a new mechanism to an existing model version we verified by extensive computer simulations (usually hundreds of runs), that within the parameter range for each parameter of the existing model no satisfying agreement between model and data could be achieved. Minimal is here to be understood as sufficient to explain the data and containing as least mechanisms as possible, whereby the building blocks of the model were chosen from those mechanisms that have already been described somewhere for any cell population. A similar iterative strategy was pursued for liver regeneration after drug induced damage predicting a previously unrecognized and subsequently validated order mechanisms [65].

We studied four different combinations of glucose and oxygen in the medium. To explain the growth kinetics, the proliferation, ECM, and cell death for the condition with high glucose and high oxygen medium concentration ([G] = 25mM, [O] = 0.28mM), the second with intermediate concentration of glucose and high oxygen concentration ([G] = 5mM, [O] = 0.28mM), we needed to assume that the cell cycle progression is possible only above a critical local production rate of ATP (= mM/h). A second necessary condition was, that the local
<figure><img src="image_14.png" /><figcaption>**Fig 10.** Waste & under-oxygenation mediating quiescence. If a cell was exposed to waste (cellular debris from cellular lyses of dead cells) or deprived from oxygen for a certain number of cell cycles, $n_{exp}^{max}$, it will become quiescent after its next division. The population kinetics of all conditions I-IV during all growth phases including the saturation as well as the spatial temporal profiles for conditions III ([G] = 25mM, [O] = 0.28mM), and II ([G] = 5mM, [O] = 0.28mM) were correctly captured. The graphs show a comparison between experimental observation and model simulation. The experimental growth curves are either shown as boxplots (box: mean, lower & upper quartiles, horizontal dashes: minimum & maximum), or as composition of mean (line) and standard deviation (error bar) for the radial profiles. The simulation curves are shown as black lines.</figcaption></figure>

doi:10.1371/journal.pcbi.1004412.g010

density of extra-cellular matrix had to be higher than a critical value (0.003). This is in accordance to literature, where dependence of cancer progression on the ECM has been shown for skin cancer [66], breast cancer [67] and NSCLC [68], where Collagen IV can regulate crucial cell signaling. If both conditions (enough ATP and ECM) were fulfilled, cells could reenter the cell cycle after a cell division. Here, cells, which were closer to the spheroid surface and thus needed less energy in order to expand, had an increasing chance to continue proliferation and not to become quiescent. Interestingly the decision whether a cell in a certain condition became
quiescent, had to be stochastic. This introduced some heterogeneity in subsets of cells in the same conditions. A deterministic scenario could not have explained the smooth transition from proliferating to quiescent zones.

The production rate of ATP depended on the local oxygen and glucose concentrations. Thereby, the ratio between both dictates to which extent a cell is in the aerobic Krebs cycle or the anaerobic lactate fermentation. Warburg stated in [69] that all cancer cells suffer from an injured respiration and thus have an exclusively anaerobic metabolism. In opposition, Zu and Guppy [70] disproved this hypothesis due to the lack of evidence and rather claimed the metabolism in cancer cells to be functional, but mainly glycolytic due to hypoxia. Here we come to a partially different conclusion: if cells are sufficiently supplied with glucose (independent from the oxygen supply), the metabolism will remain glycolytic (90%), and only if the glucose supply is getting short, the metabolism will favor the aerobic Krebs cycle (see Fig 8). Besides lactate acidity (> 20mM), the depletion of carbon sources to maintain a critical ATP production (= 900mM/h) and not hypoxia were the main reasons of death. The latter was also recently suggested by Kasinskas et al. [71], while, in contrast to our assumptions, they excluded lactate as source of acidity and instead assumed it to be an important secondary metabolic resource. However, either growth adverse or death promoting effects were described for high lactate concentrations [59]. So here further clarification of the dominating role of lactate would be necessary. The functional forms of the oxygen and glucose consumption rates were inferred from experimental findings of Freyer, Sutherland and co-workers in EMT6/Ro cells. The lactate and ATP production rates were then directly derived from those rates by the single assumption that cells transform the consumed glucose in an optimal way with respect to ATP output. For wide ranges of glucose and oxygen concentrations the ATP production rates remain stable between 80...130 × 10⁻¹⁷ mol/cell/s or 1000...1700mM/h respectively, assuming a reference cell volume of 2700µm. In literature values can be found between 4.6...15.3 × 10⁻¹⁷ mol/cell/s ([72-75]). The difference could be either due to differences in energy needs between different cell types, or to the model simplification that glucose in our model is exclusively used for metabolism.

Interestingly and importantly, the model, despite only having been calibrated with two of the four growth conditions, were subsequently able to correctly and quantitatively predict the growth phase of the other two growth conditions ([G] = 1mM, [O] = 0.28mM and [G] = 25mM, [O] = 0.07mM, respectively). This indicates that the model did capture the functionalities necessary to explain the data for different glucose and oxygen conditions. To further permit independent validation of our model, we performed additional simulations for other glucose and oxygen medium concentrations (Fig S12 in S1 Document).

However, all growth curves showed saturation and partially even shrinkage after some time. The saturation phase could be largely captured by adding the potential effect of a waste produce being released in the extracellular space from cells undergoing lysis. Shrinkage could be added if dying cells at the border detach and enter the growth medium; however, we did not consider this process, as it was not observed in the experiments (for example, for A549 cells, another NSLC cell line, a massive detachment of cells from the spheroid could be observed in the experiments). Interestingly, model simulations with a lysis rate of 0.35/h, a typical value in-vivo, turned out to be incompatible with the in-vitro data. A lysis rate of a few hours as observed in-vivo would lead to a very fast removal of dying cells and thus almost no dead cells in the tumor center, in sharp contrast to the in-vitro experiments. We obtained a much smaller value of about 0.01/h by comparison of model simulation results and the spatial cell death and proliferation profiles i.e., only such a small lysis rate permits the occurrence of a “necrotic core” as observed in the in-vitro experiments. For such a low lysis rate we found that the apoptosis—if present—would need to be very slow, as it affects also cells in the viable rim in order
to agree with the experimental observation of only very little dead cells in the viable rim. For this reason, apoptosis could be neglected in explaining the experimental results in this paper. The small value of the lysis rate, even though surprising on a first view, may be explained by noticing that stromal cells (such as e.g. macrophages) digesting dead cells are not present in-vitro. Hence lysis might be expected to be slower in-vitro than in-vivo.

Contact inhibition seems to be a crucial element. Suppressing contact inhibition with varying combinations of the other mechanisms in each case leads to complete failure of match between data and model simulations (see Fig S6 in S1 Document, where the parameters of model 4 has been used). This observation supports the view expressed previously in the paper that a mechanical growth inhibition plays an important role in multicellular spheroids.

We moreover tested the possibility that cells may actively migrate towards the necrotic zone by necrotaxis (Figs S9, S10 in S1 Document). As to keep a sufficiently large necrotic core as experimentally observed the lysis rate had to be small, significant migration could not be observed. On the other hand, if the lysis rate was chosen large, then significant migration of cells could be observed but the necrotic core was too small, as cells in the center were too quickly eliminated by lysis. In the latter case, the necrotic core with increasing migration rate became smaller (Figs S9, S10 in S1 Document). We concluded that migration driven by morphogens towards the central necrosis in SK-MES-1 cells is small.

Interestingly the final model emerging from this stepwise, image-guided inference strategy closely resembles the hypothesis on growth control of MCTS by growth promoters (GP), growth inhibitors (GI), viability promoters (VP) and inhibitors (VI) (Fig 11).

In order to permit validation of our model, we simulated a number of predictions (Fig S13 in S1 Document). We predicted the spatial temporal growth dynamics for [G] = 1mM, [G] = 3mM, O₂ = [0.28mM], and [G] = 25mM, [O] = 0.07mM. In this context we would like to remind that our model was able to predict the growth kinetics (L(t)) for [G] = 1mM, [O₂] = 0.28mM, and [G] = 25mM, [O₂] = 0.08mM correctly.

<figure><img src="image_15.png" /><figcaption>**Fig 11.** Spheroids show a radial organization of cell phenotypes: proliferating (green), quiescent (blue) and dying cells (black). The thickness of those layers is controlled by different factors as growth promoters (GP), growth inhibitors (GI), viability promoters (VP) and inhibitors (VI). The different control molecular components can largely be mapped on the control molecules emerging in our final model.</figcaption></figure>

doi:10.1371/journal.pcbi.1004412.g011
In order to quantify the goodness of the fits shown in Figs 7, 9 and 10 we calculated the log-likelihoods for subsets of curves (see Table 3A) as well as the ensemble of all curves (see Table 3B). We assumed that the measurement error is additive, normally distributed as well as independent and identically distributed (i.i.d.). Accordingly, the likelihood of the measured data mean μ given the parameter θ and the corresponding model prediction x(θ) is given by
$$L(\theta) = \prod_i \frac{1}{\sqrt{2\pi\sigma_i^2}} \exp\left(\frac{(x(\theta)_i - \mu_i)^2}{2\sigma_i^2}\right),$$
where the index i runs over the data points. The uncertainty of the data points is determined by the standard deviation, σ. Accordingly, points with large uncertainties σ are weighted less. Despite some deviations in single profiles, the final model considering ATP, lactate and waste is found to be the most likely to explain the experimental data. On the other hand, an increase of the likelihood correlates with the increasing number of model parameters and the risk of over-fitting. Especially model 2 and 3 have a small relative difference in log-likelihood. As a measure for model quality accounting for the number of parameters we used the Akaike information criterion (AIC) (see Table 3D), which also confirms the final model to be the best choice. We note that the number of parameters had no influence on the ordering of the models as the absolute differences of their log-likelihoods ΔlnL is many orders of magnitude larger than the difference in number of parameters Δk (see Table 3B and 3C).

To conclude, we would like to stress the key message demonstrated in this paper: quantitative comparison of spatial profiles observed from cells states in different experimental conditions and time-points, generates information so rich that one may infer even molecular control mechanisms and parameters of spatio-temporal growth and death patterns. Hence, careful imaging, image processing and image analysis may serve as an important source of information to infer mechanistic knowledge on tissue growth and organization processes. Such an approach would gain to be more fully explored.

Our model is hybrid. It integrates as separate components the cell and molecules, and as functional components a mechanical form of contact inhibition, a metabolic component comprising oxygen, glucose, ATP, lactate, and waste, several of the molecules acting as morphogens. It would be interesting to see in how far the same model can capture the growth behavior of other cell lines and of other cell types.

<table><caption>**Table 3. Model comparison.** The table shows the log likelihoods (A, B) and the Akaike information criterion (D) comparing three models to experimental data. The growth curves of all nutrient conditions and the radial profiles of proliferation, death and ECM at a specific time were compared individually (A) as well as a whole (B, D). We also calculated the values according to the Bayesian information criterion (BIC), which led to the same ordering (not shown). Notice that the penalty terms in the AIC (2k) and BIC (k log(n), n: number of data points) play no role in the evaluation as the difference in the penalty terms are much (4–5 orders of magnitude) smaller than ln L. (Δ AIC = AIC(Model row-index)—AIC(Model column-index).)</caption><thead><tr><th></th><th>Log-likelihood<br>In L</th><th>(Model 2)<br>ATP</th><th>(Model 3)<br>ATP & Lactate</th><th>(Model 4)<br>ATP, Lactate & Waste</th><th>#Data Points<br>n</th></tr></thead><tbody><tr><td>A</td><td>growth curves (all conditions)</td><td>-8.61e+02</td><td>-7.42e+02</td><td>-1.72e+02</td><td>60</td></tr><tr><td></td><td>proliferation (conditions II & III, day 17)</td><td>-1.83e+05</td><td>-1.60e+05</td><td>-1.83e+05</td><td>460</td></tr><tr><td></td><td>death (conditions II & III, day 17)</td><td>-2.69e+03</td><td>-5.15e+03</td><td>-6.16e+03</td><td>532</td></tr><tr><td></td><td>ECM (conditions II & III, day 17)</td><td>-1.23e+06</td><td>-1.23e+06</td><td>-8.66e+05</td><td>524</td></tr><tr><td></td><td>proliferation (conditions II & III, day 24)</td><td>-1.97e+02</td><td>-1.40e+02</td><td>1.18e+02</td><td>532</td></tr><tr><td></td><td>death (conditions II & III, day 24)</td><td>-7.47e+02</td><td>2.16e+01</td><td>-3.11e+02</td><td>616</td></tr><tr><td></td><td>ECM (conditions II & III, day 24)</td><td>-1.82e+01</td><td>3.08e+01</td><td>3.26e+02</td><td>609</td></tr><tr><td>B</td><td>total</td><td>-1.42e+06</td><td>-1.39e+06</td><td>-1.06e+06</td><td>3333</td></tr><tr><td>C</td><td>#model parameters k</td><td>reference</td><td>+2</td><td>+9</td><td></td></tr><tr><td>D</td><td>(Model 2)</td><td>0.00E+00</td><td>-5.00E+04</td><td>-7.30E+05</td><td></td></tr><tr><td></td><td>ΔAIC = 2Δk - 2Δln L (Model 3)</td><td>5.00E+04</td><td>0.00E+00</td><td>-6.80E+05</td><td></td></tr><tr><td></td><td>(Model 4)</td><td>7.30E+05</td><td>6.80E+05</td><td>0.00E+00</td><td></td></tr></tbody></table>

doi:10.1371/journal.pcbi.1004412.t003
We think that the possible imaging techniques and image analysis software in combination with modeling could permit a screening of growth dynamics and subsequent quantitative classification of multicellular spheroids. For example, EMT6/Ro cells ([14]) show a very similar but not equal growth phenotype as SK-MES-1 cells: (1) detachment of cells is rare (as opposed to, for example, A549 cells, that reveal significant detachment in-vitro (and in-vivo)), (2) under sufficient oxygen supply, the growth of the outer spheroid diameter remains unaffected (or almost unaffected) by glucose, while (3) reduction of oxygen from 0.28mM to 0.07mM reduces growth dramatically in SK-MES-1 cells even if glucose medium stays high, while in EMT6/Ro cells no reduction is observed as long as nutrient medium concentration stays high: this demarcates a difference between the EMT6/Ro and SK-MES-1 cells. On the other hand, glucose affects the size of the necrotic core. Reduction from 16.5mM glucose to 0.8mM glucose medium concentration in EMT6/Ro cells (at 0.28mM oxygen) increases the necrotic core (as one can infer from comparing the cell count with the diameter, see [34]), which can also be observed in SK-MES-1 cells if glucose is reduced from 25mM to 5mM (at 0.28mM O₂). However, in SK-MES-1 cells the necrotic core for richer glucose ([G] = 25mM vs. 5mM) occurs later but at about 24 days is about the same size.

Models can provide a quantitative framework to test how far such differences can be attributed to parameters with the same model, or whether “another” model needs to be used by adding or dropping mechanisms. For example in the first case, can the same model be used to capture a wide range of cell lines with regard to their MCS growth behavior by only adjusting its parameters—indicating only quantitative changes, or, in the 2nd case, does one need to implement mechanisms for one cell type that are not observed within the physiological range of parameters for another cell type? Given how much multicellular spheroids are still in use as biological model system, we think it would be of fundamental interest to do such an analysis as a community effort, even though this might be considered as on the first view as a “step-back” as the growth dynamics of multi-cellular spheroids could have been measured 20–30 years ago. Modern technology could largely permit automated analysis if pipelines were constructed for that purpose, hence avoiding the largely manual and tedious analysis applied for the work in this paper. Adding more and more cell lines would permit to refine the model one starts with, and zoom into the so far still highly simplified representation of mechanisms without the threat of having a far to large number of fit parameters that cannot be controlled. In this way, identification of necessary model components and adjustment of parameters linking the components could be achieved.

# Materials and Methods

## Cell Culture

NSCLC cell line SK-MES-1 used in this study was obtained from ATCC (Manassas, VA, USA) and cultivated in a humidity controlled incubator at 37 C and 5% CO₂ in 150cm² tissue culture dishes (TPP) in DMEM (Dulbecco's modified Eagle's medium, LONZA, Verviers, Belgium) supplemented with 10% FCS (fetal calf serum, Southern America, GIBCO, Germany) and 1% Penicillin/Streptomycin (Biochrom AG, Berlin, Germany). Cells were used between passages 10 and 30 and passaged at a split ratio of 1:4 to 1:6. Cultures were routinely tested for mycoplasm contamination as described by Stacey and Doyle 1997 and always found to be negative. Additional medium for the test cultures was DMEM w/o Glucose (GIBCO, Germany) supplemented with 10%FCS and 1mM Glucose (Carl Roth GmbH, Germany), and DMEM with 1.0 g/L glucose w/o L-Glutamine (LONZA, Verviers, Belgium) supplemented with 10% FCS and 25mM L-Glutamin (SIGMA, Germany). Additionally, cells were kept in a humidity controlled incubator at 37 C and 5% CO₂ and either normal atmospheric 20% O₂ (corresponding to
0.28mM) or 5% O₂ (corresponding to 0.07mM). This choice permits comparison to classical work in literature (e.g. [12]) and takes into account that lung is rich in oxygen. NSLC cells originate from lung epithelium having at least partially direct contact to the inhaled air so are at least initially not limited to the blood oxygen level.

**Solutions.** 2% Methocel solution was prepared by stirring 6g Methylcellulose (SIGMA, Germany) in 250ml propagation medium at 60 C for 20min. Then a further 250ml of medium were added and stirred at 4 C over night. The solution was aliquotted in 50ml Falcon tubes, centrifuged two times for 99min at 4000 rpm to concentrate long Methylcellulose-Fibers at the bottom of the tube. Tubes were then stored upright at 4°C.

2% Agar-Solution was prepared by solving 5g bacterial grade Agar (GIBCO, Germany) in 250ml of H2O. This solution was autoclaved and then kept at 60 C until 24 well-Plates were coated.

**Spheroid generation and cultivation.** For the generation of spheroids with defined size and cell number, a hanging drop assay was employed. Here cells were first trypsinized, counted, then centrifuged and resuspended in propagation medium with 20% Methocel-Solution. Drops of 20µl of this cell suspension were pipetted on the lid of a 150cm²-culture dish, which was subsequently carefully inverted back onto the dish. In this approach, all suspended cells in the resulting hanging drop contribute to the formation of a single spheroid. After 48h in the hanging drop the spheroids were transferred on 24 well plates, by carefully pipetting with a cut 200ml pipette tip, with one spheroid in one ml of respective medium per well. The wells were pre-coated with 250µl of 2% agar each to prevent attachment of spheroids. Medium was changed once a week. Four different glucose/oxygen conditions were employed (see Table 1). Spheroids were cultivated over a period of over 48 days.

## Data Acquisition

**Spheroid growth curves.** Growth of spheroids was monitored by acquisition of bright field images through an Olympus IX-70 microscope fitted with an AxioCam Erc5s camera (Zeiss, Germany) twice a week. The projected area of spheroids was determined using ImageJ, software version 1.43u. Mean areas and standard deviations were calculated using Microsoft Excel. At least 4 spheroids were evaluated per time point and condition.

**Cryosectioning and immunofluorescence staining.** Spheroids were taken out at specified time points (see Table 4), embedded in TissueTek (SAKURA Finetek, Netherlands) cryo-medium, frozen over liquid nitrogen and processed for cryostat sectioning. Cryosections of 6–8 µm thickness were mounted on slides, air dried and then fixed with 4% PFA for 20 min at RT, washed in PBS for 30 min at RT, permeabilized with ice cold 0.1% TritonX-100 / 0.1% sodium citrate for 2min on ice. To stain for dead cells we employed the In Situ Cell Death Detection Kit (No. 12 156 792 910 Roche Applied Science) according to the manufacturer’s protocol. Staining for proliferating cells and ECM component Collagen IV was performed with anti human Ki67, mouse monoclonal (No. M7240 DakoCytomation, Glostrup, Denmark) and anti Collagen IV, rabbit polyclonal (No. 10760, Progen, Heidelberg, Germany) antibody in 12%

<table><caption>**Table 4. Timepoints of cryosectioning.**</caption><thead><tr><td>Name</td><td>Time</td></tr></thead><tbody><tr><td>T3</td><td>17days</td></tr><tr><td>T4</td><td>24days</td></tr><tr><td>T5</td><td>34days</td></tr><tr><td>T6</td><td>46days</td></tr></tbody></table>

doi:10.1371/journal.pcbi.1004412.t004
BSA respectively. Secondary antibodies used were anti mouse, donkey, Cy3 (No. 715-166-151, Dianova, Hamburg, Germany) and anti rabbit, goat, Alexa 488 (No. A-11034 Molecular Probes USA/NL) in 12% BSA. In both cases, nuclei were stained with HOECHST (Bisbenzimide H 33258 SERVA, Heidelberg, Germany) in 12%BSA. Sections were examined and photographed using a Leica DM RBE (Leica, Germany) microscope fitted with epifluorescence optics.

## Quantitative Image Analysis

The images acquired (see above) are raw data. They consist of a set of pixels with position $(x, y) \in \mathbb{N}^2$ and color intensities for the three color channels red, green and blue defined by

$$ 
\begin{aligned}
I^{\text{channel}}[x, y] &\in [0, 1), \\
I^{\text{channel}}: \mathbb{N}^2 &\to \mathbb{R},
\end{aligned} 
\qquad (12) 
$$

where channel ∈ {red, green, blue}. In the following, the tools used to preprocess the raw images (e.g. to reduce noise) and to identify or segment objects (e.g. cell nuclei, spheroid border) will be introduced. Then the preprocessed images have been analyzed quantitatively as described subsequently.

**Nuclei segmentation.** The cell nuclei were stained with HOECHST staining (blue color channel, see Fig 2). In order to segment the single nuclei, the images were smoothed by a series of four median filters with kernel size 3 × 3 (see below) as a preprocessing step. A larger kernel would smooth the image too much. In the following, the watershed algorithm(see below) was applied to the inverted blue color channel, 1 − $I[x, y] (I[x, y] \in [0, 1))$ for pixels above a certain threshold $I[x, y]^{blue} > I^{HOECHST}$ (see Table 4). Introduction of a threshold avoids the segmentation of the whole image, focusing only on regions clearly stained by HOECHST and thus associated with the cellular nuclei.

**Segmentation (watershed).** The watershed algorithm is a segmentation algorithm separating an image into different regions using the concept of watersheds. It was first proposed in ref. [76], and later enhanced by a fast algorithm presented in ref. [77]. The idea of the algorithm is to treat the images as landscapes, where color intensities $I[x, y]$ correspond to the amplitude. The landscape representation of an example image is shown in Fig 4. The local minima of the landscape are identified with catchment basins where water would accumulate. To detect the pixels belonging to the different basins, the landscape is successively “flooded” and pixels of the current altitude are associated to the basin of their neighboring pixels (already associated to a basin). When two basins touch each other a watershed line between them is created and expansion is stopped. The algorithm either stops at a chosen maximum altitude (color intensity), or when none of the basins can expand anymore.

**Noise-reduction (median filter).** As each local minimum is giving rise to one basin, it is important to smooth the images thereby reducing noise in the image before performing the segmentation in order to avoid an over-segmentation into many fragmented pieces. For this purpose a number of different linear (e.g. mean filter, Gaussian blur) and non-linear filters (median) were applied. We used a sequence of median filters known to result in an edge-preserving smoothing, i.e. while removing noise it conserves the original features of the image (e.g. nuclei shapes) in contrast to, for example, Gaussian blur. The impulse response $I'[x, y]$ for each input pixel $(x, y)$ is the median of neighboring entries within a certain ‘window’ of size $n \times n$:

$$ I'[x, y] = \text{median}_{|x-x'| \le a, |y-y'| \le a} \{I[x', y']\}, \qquad (13) $$

where $a = (n - 1)/2$.
<figure><img src="image_16.png" /><figcaption>**Fig 12. Cell density and cell size estimation.** Left: The Delaunay triangulation of all segmented nuclei serves to estimate the cell sizes via its dual, the Voronoi diagram. Right: Average cell diameter as a function of distance to the spheroid border. The black curve is the average profile of six images (condition III, T3) with bin size 1µm and the red curve is the gliding average with window size 10µm.</figcaption></figure>

doi:10.1371/journal.pcbi.1004412.g012

**Spheroid lumen and border.** The basis for estimating the spheroids borders was the cell nuclei (see above). By applying a series of dilation filters a coarse approximation of the cell shapes around the nuclei was achieved. A dilation can be described as the following:

$$I'[x, y] = \max_{|x-x'| \le a, |y-y'| \le a} \{I[x', y']\}.$$

Finally, the empty spaces (necrotic areas) in the center of the spheroids were removed by hole-filling algorithms (e.g. by flood filling). Fig 4 shows the estimated spheroid lumen around the cell nuclei in green color for an example picture. The spheroid border is assumed to be the interface between the green (lumen) and the black (background) area(s).

**Radial profiles, binning and averaging.** As the spatial arrangement of cell phenotypes (proliferating, quiescent and dead) and ECM density could be assumed to be approximately concentric, the following sections focus on how to extract the spatial statistical profiles of the quantities stained. As Fig 4 shows, the spheroids were not completely spherical, but typically slightly elongated with irregular surfaces. To account for this, all statistics have been performed as a function of distance to the spheroid border, defined by $\Delta L$, as opposed to the distance to the center of mass, which would be the method of choice in case the MCTSs would have been perfectly spherical. For estimation of the radial profiles binning was utilized. I.e. $\Delta L$ was divided into small intervals $\Delta L_i = [h(i), h(i + 1))$, so-called bins, of $h = 1\mu m$ length (image resolution: pixel = 0.98$\mu m^2$). Then each nucleus that contains at least one pixel, whose distance to the closest border pixel is $\Delta L \in \Delta L_i$, enters into the statistics of bin $i$. The average curves were generated from the radial profiles of the individual images by calculating the average values of the corresponding bins with the same index.

**Nuclei density and cell diameter.** The local cell density can be estimated from the number of nuclei per area. Its inverse is the area per nucleus. The nuclei served as construction points for a Voronoi diagram. In a Voronoi diagram all points closer to one construction point than to any other construction point belongs to one Voronoi cell. Voronoi cells are polygonal in shape. Accordingly, the Voronoi cell belonging to one nucleus was identified with the 2D
<figure><img src="image_17.png" /><figcaption>Fig 13. Automatic versus manual detection of Ki67 positive nuclei.</figcaption></figure>

doi:10.1371/journal.pcbi.1004412.g013

A good agreement (small classification error) between automated and manual Ki67 detection was reached for $0.35 \le I^{Ki67} \le 0.5$ and $\phi^{Ki67} \le 0.2$. The parameter values used in the following are denoted in Table 5 and the corresponding comparison of manually and automatically detected Ki67-positive nuclei is shown in Fig 13.

## Cellular Automaton

The cellular automaton model, which extends on previous work ([43, 47]), is defined by a set of rules:

**Cell.** One biological cell can occupy either one or two sites on a Voronoi diagram (Fig 5, Fig S1 in S1 Document). The term Voronoi lattice site is used from now on to distinguish the Voronoi cell from a biological cell. At the beginning of the cell cycle, i.e. after division, it always occupies only one Voronoi lattice site. Later in the cell cycle (specified below) it grows by occupying a neighboring lattice site as well.

**Cell cycle progression/replication.** The cell cycle is subdivided into $m_d$ intervals (Fig 5). At the beginning of the cell cycle the cell internal counter *M* is set to *M* = 0. If the nutrient conditions are permissive for proliferation according to the conditions detailed in the main text
<figure><img src="image_18.png" /><figcaption>**Fig 14.** Sensitivity analysis for threshold parameters $I^{K167}$ and $\phi^{K167}$. The plots show the sensitivity, specificity and classification error comparing the manual detection with the automated detection for a combination of {$I^{K167}$, $\phi^{K167}$}.</figcaption></figure>

doi:10.1371/journal.pcbi.1004412.g014

and in the supporting information (S1 Document) then a cell enters the cell cycle. This was mimicked by the transition $M = 0 \rightarrow M = 1$. A proliferating cell successively increased its internal counter starting from $M = 0$ with a rate $m_d k_{div}$ by one until $M = m_d$. Here $k_{div} = \tau^{-1}$, where $\tau$ is the expectation value of the cell cycle time, and depends on the environmental conditions

<table><caption>**Table 5.** Image processing parameters.</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>μ<sup>HOECHST</sup></td><td>0.15</td></tr><tr><td>μ<sup>TUNEL</sup></td><td>0.45</td></tr><tr><td>φ<sup>TUNEL</sup></td><td>0.20</td></tr><tr><td>μ<sup>K167</sup></td><td>0.35</td></tr><tr><td>φ<sup>K167</sup></td><td>0.05</td></tr></tbody></table>

doi:10.1371/journal.pcbi.1004412.t005
(see Eq 11). Note that the transitions are Poisson processes; the emerging total distribution of cell cycle time is an Erlang distribution [43].

**A) Cell growth.** At $M = m_g \le m_d$ a cell grows by occupying a neighbor Voronoi lattice site in addition to the one it occupies already. Consequently, for the choice $m_g < m_d$ the doubling of the mass took place somewhere in the cell cycle but not in the moment of cell division as it is the case in the classical Eden model [78] or more recent work [43]. If all neighboring lattice sites are already occupied, then the cell first liberates one of them by pushing adjacent cells toward the closest free lattice site and then occupies the original site and the liberated neighbor site (for details see supplementary material).

**B) Cell division.** When $M = m_d$ the cell splits into two daughter cells with each of the daughter cells occupying a single Voronoi lattice site. After the division each cell decides whether to reenter the cell cycle or to become (permanently) quiescent. The probability to reenter the cell cycle, $p_{re}$, depends on the environmental conditions (see Eq 4).

**Cell death.** Under certain conditions specified in the main body of the text a cell dies with the rate $k_{nec}$ (see Eq 9).

**Cell lysis.** Dying cells are lysed with the rate $k_{lys}$. Consequently the corresponding Voronoi lattice sites are set free.

**Cell migration.** Cells are able to move either by being pushed if another cell grows, or by active migration. Active migration is mimicked by hopping to a neighboring Voronoi site. The choice of the hopping rate depends on the considered biological process:

**A) Random walk/free diffusion.** In the case of unbiased cell movement the hopping rate has been chosen to $\lambda_i = \lambda/n$ for all n neighboring sites $i = 1...n$. This corresponds to a diffusion process such that $\lambda = 6D/l^2$, with D the diffusion coefficient and l the average distance between two adjacent lattice sites.

**B) Biased random walk.** The dynamics of the multi-cellular configuration is determined by a master equation for the multivariate probability to find a certain multi-cellular configuration, X where $X = \{x_1, x_2, ...\}$ denotes the state vector of the multicellular configuration. One way is to enumerate all lattice sites and denote by $x_k$ the state vector of the cell localized at lattice site k. If this lattice site is empty, the state is zero. The dynamics is then formalized by

$$ \frac{\partial p(X, t)}{\partial t} = \sum_{x'} e^{\frac{-\Delta E_{X \to X'}}{F_T}} \cdot (\tau_{X \to X})^{-1} \cdot p(X', t) - e^{\frac{-\Delta E_{X \to X'}}{F_T}} \cdot (\tau_{X \to X'})^{-1} \cdot p(X, t). $$

In the case of cell-cell adhesion $\Delta E_{X \to X'}$ represents the energy difference resulting from a change of cell-cell contacts when moving from the current (configuration X) to the new lattice site (configuration X'). $p(X, t)$ is the probability of the system being in configuration X at time t. $\tau_{Y \to Z}$ denotes the time between two successive hops from Y to Z. $F_T$ is a characteristic fluctuation energy [1]. An alternative explained in the supplementary information is to enumerate all cells, and for each cell the processes it can perform. A process is then chosen with a probability that corresponds to its relative weight, calculated as the rate for this process divided by the sum of the rates for all other processes (for algorithmic details, see supplementary information).

## Metabolism

Glucose and oxygen are among the main metabolites of most biological cells. In normal tissue they are provided mainly by the vascularization. Tumor spheroids as avascular tumors (in-vivo) are mainly fed by nutrients diffusing into the interior from the border. Above a certain
size they display regions lacking glucose (hypo-nutrition) and/or oxygen (hypoxia). This occurs if the nutrients entering the tumor via its borders are consumed completely before reaching the center.

We modeled the diffusion and uptake of glucose and oxygen by a reaction-diffusion equation:

$$ \partial_t u = \nabla \cdot (D_u(\sigma)\nabla u) + r_u(\sigma, u), \qquad (19) $$

where $D_u$ is the molecule diffusion coefficient and $r_u(\sigma, u)$ denotes the reaction term, $u \in \{G, O\}$. The diffusion coefficient was chosen differently in the nutrient medium than inside the tumor spheroid. The reaction term mimicked the cells consumption with the choice

$$ r_u(\sigma, u) = -q_u(u)\sigma = -q_u(u)\sum_k \delta(\mathbf{r} - \mathbf{r}_k). \qquad (20) $$

**Glucose and oxygen consumption rates.** In experiments on EMT6/Ro cells for varying medium concentrations the cell consumption rates for each of the two molecules were observed to depend on each other ([12, 79]). Casciari et al. [80] proposed an extended Michaelis-Menten-like consumption rate for glucose,

$$ q_G = V_G^{max} \frac{[G]}{[G] + k_G^G}, \qquad (21) $$

and oxygen,

$$ q_{O_2} = V_{O_2}^{max} \frac{[O_2]}{[O_2] + k_{O_2}^{O_2}}. \qquad (22) $$

<figure><img src="image_19.png" /><figcaption>**Fig 15.** Oxygen (right) and Glucose (left) consumption. The graphs show a comparison between experimental measurements of glucose and oxygen consumption rates of EMT6/Ro cells [80] (dots) and predictions of Eqs 21 and 22 (grid). The rates are functions of the local extra-cellular glucose and oxygen concentrations. (The original values were rescaled from mol/cell/s to mM/h (= mol/m³/h) for an average cell volume of 2700µm: 1mol/cell/s = 13.3 x 10¹⁷ mM/h.)</figcaption></figure>

doi:10.1371/journal.pcbi.1004412.g015
By direct comparison to their data we found the following cross-linking terms to fit well the experimentally observed consumption rates (Fig 15):

$$V_G^{\max} = q_G^{\max} \left( 1 - \left( 1 - \frac{q_G^{\min}}{q_G^{\max}} \right) \frac{[O_2]}{[O_2] + k_G^{O_2}} \right), \qquad (23)$$

$$V_{O_2}^{\max} = q_{O_2}^{\max} \left( 1 - \left( 1 - \frac{q_{O_2}^{\min}}{q_{O_2}^{\max}} \right) \frac{[G]}{[G] + k_{O_2}^{G}} \right).$$

Within this paper we assumed the same functional form and parameter values for the consumption rates of the SK-MES-1 cells as no equivalent measurements as those for the glucose and oxygen consumption rates of EMT6/Ro cells were available for SK-MES-1 cells. The parameters of Eqs 21, 22 and 23 were inferred from [12] (see Fig 15 and Table 6).

**ATP as energy currency of the cell.** As glucose or oxygen-dependent thresholds for cell cycle progression or cell death led to results disagreeing with the experimental observations, in a next step we investigated the potential impact of glucose and oxygen on cell division and death. Besides being utilized as cell material, glucose is mainly used to gain energy for cell maintenance and reproduction. This happens by oxidation of glucose, which under standard conditions reads [54]:

$$C_6H_{12}O_6 + 6O_2 \rightarrow 6CO_2 + 6H_2O.$$

The freed energy ($\Delta G = -2880kJ/mol$) is thereby stored in ATP. We below denote in a highly simplified form by summarizing reactions the contributions of glucose and oxygen to the ATP production. For one mol of glucose, within Glycolysis,

$$\text{Glucose} + 2\text{NAD}^{+} \xrightarrow{k_1} 2\text{Pyruvate} + 2\text{ATP} + 2\text{NADH}. \qquad (24)$$

Under anaerobic conditions (lack of oxygen), pyruvate is metabolized into lactate in order to convert NADH back to NAD$^+$ (see Eq 29).

Otherwise, if sufficient oxygen is available, the two moles of NADH are turned into 4-6 ATP molecules by oxidative phosphorylation, depending on the shuttle system transporting the electrons of the NADH through the mitochondrial membrane from the cytosol into the mitochondria (the inner mitochondrial membrane is impermeable for NADH or NAD$^+$). For glycerin-3-phosphat-shuttle it is 4 ATP. While glycolysis occurs in the cytosol, the citric acid cycle occurs in the mitochondrial matrix, and oxidative phosphorylation at the inner mitochondrial membrane. In the presence of sufficient oxygen, under aerobic conditions, one mol

<table><caption>**Table 6. Model parameters for glucose and oxygen consumption.**</caption><thead><tr><th>Parameter</th><th>Unit</th><th>Value</th></tr></thead><tbody><tr><td>k<sub>G</sub><sup>G</sup></td><td>mM</td><td>0.068</td></tr><tr><td>k<sub>G</sub><sup>O<sub>2</sub></sup></td><td>mM</td><td>0.031</td></tr><tr><td>q<sub>G</sub><sup>min</sup></td><td>mM/h</td><td>186.67</td></tr><tr><td>q<sub>G</sub><sup>max</sup></td><td>mM/h</td><td>706.67</td></tr><tr><td>k<sub>O<sub>2</sub></sub><sup>O<sub>2</sub></sup></td><td>mM</td><td>0.031</td></tr><tr><td>k<sub>O<sub>2</sub></sub><sup>G</sup></td><td>mM</td><td>0.100</td></tr><tr><td>q<sub>O<sub>2</sub></sub><sup>min</sup></td><td>mM/h</td><td>120.00</td></tr><tr><td>q<sub>O<sub>2</sub></sub><sup>max</sup></td><td>mM/h</td><td>306.66</td></tr></tbody></table>

doi:10.1371/journal.pcbi.1004412.t006
of pyruvate is transformed to one mol of acetyl-CoA, and introduced into the citric acid cycle whereby one mol of NADH is generated. In the citric acid cycle, from one mol of acetyl-CoA, 3 moles of NADH, one mole of FADH₂, and one GTP are generated. By oxidative phosphorylation, one mol of NADH yields 3 moles of ATP, and one mol of FADH₂ yield 2 moles of ATP. For each mol of NADH or FADH₂ entering oxidative phosphorylation, 1/2 mol of O₂ is consumed. Summarizing, in the presence of sufficient oxygen, for one mol pyruvate, entering the citric acid cycle, 4 moles of NADH yielding 12 moles of ATP, 1 mol of FADH₂ yielding 2 moles of ATP, and 1 mol GTP energetically corresponding to one mol of ATP are produced. As NADH from glycolysis (see Eq 24) is only turned to ATP in presence of sufficient oxygen, we add the emerging 2–3 ATP per NADH to the r.h.s. and the NADH to the l.h.s. of the aerobe pathway reaction:

$$ \text{Pyruvate} + \text{NADH} + 3\text{O}_2 \xrightarrow{k_2} 17(18)\text{ATP} + \text{NAD}^+. $$

As mentioned above, the reactions take partially place in the cytosol, partially in the mitochondria, and the ATP yield from NADH depends on the shuttle system used which is why the numbers of ATP yield can slightly differ. We here continue our calculation with 17 moles ATP per mol pyruvate.

We assume that the changes in intra-cellular glucose and oxygen concentrations mainly depend on the metabolic rates, $k_1$ and $k_2$, and the uptake rates from the extra-cellular space, $q_G$ and $q_{O_2}$,

$$ \frac{d[G]}{dt} = -k_1[G] + q_G, \qquad (25) $$

$$ \frac{d[\text{O}_2]}{dt} = -3k_2[\text{Pyruvate}][\text{NADH}][\text{O}_2]^3 + q_{\text{O}_2}. \qquad (26) $$

The changes in ATP concentration can be written as following:

$$ \frac{d[\text{ATP}]}{dt} = 2k_1[G] + 17k_2[\text{Pyruvate}][\text{NADH}][\text{O}_2]^3. \qquad (27) $$

Assuming the equilibrium conditions of Eqs 25 and 26 are reached very fast we can reformulate Eq 27, and in case of $[G] > [G]^{crit}$ obtain

$$ \frac{d[\text{ATP}]}{dt} = 2q_G + \frac{17}{3}q_{O_2} = p_{\text{ATP}}. \qquad (28) $$

In case $[G] \le [G]^{crit}$, the concentration of pyruvate becomes zero, hence no ATP is produced. To avoid that cells completely deprived of glucose are still capable of producing ATP, we assume the critical amount of glucose to be $[G]^{crit} = 0.01mM$.

Fig 8 illustrates how cells regulate their consumption of glucose and oxygen under very different circumstances, hyponutrition (low glucose) and hypoxia (low oxygen), in such a way that a yield of ATP is maintained at values between 1000 and 1700 ×10⁻¹⁷ mM/h per cell. This observation suggests that continuously dividing EMT6/Ro cells need at least 1000mM/h. On the other hand Eq 28 implies that EMT6/Ro cells metabolize only 10% of the consumed glucose in the Krebs-cycle under sufficient glucose and oxygen supply (see Fig 8). The other 90% are either fermented to lactate or used to build up other cell components (e.g. amino acids, nucleotides) nevertheless important for cell reproduction. Jiang et al. [22] indicated a much lower fraction (75%) of anaerobically metabolized glucose.
**Lactate—the side product of anaerobe metabolism.** Under anaerobic conditions the cell accumulates pyruvate molecules generated in glycolysis as these cannot be further processed in the citric acid cycle (or Krebbs cycle) without the presence of oxygen. On the other hand, NADH cannot be oxidized back to NAD⁺, which is needed for the glycolysis to continue. To over come this shortcoming of NAD⁺, pyruvate is reduced to lactate by incorporating a proton of NADH.

$$ \text{Pyruvate} + \text{NADH} \xrightarrow{k_3} \text{Lactate} + \text{NAD}^+. \qquad (29) $$

This process is called the lactic acid fermentation.

In the model we assumed lactate to be a direct side product of glucose metabolized anaerobically i.e., in absence of oxygen. For each molecule of glucose entering the anaerobic lactate fermentation the cell will produce 2 molecules of lactate. Assuming the lactate production rate to be $p_L \approx 2q_G$ (where $q_G$ corresponds to anaerobic rate of glucose metabolisation) and lactate to leak out of cells and diffuse in the tissue, results in the following formulation of the lactate dynamics

$$ \frac{\partial[L]}{\partial t} = \nabla \cdot (D_L(\sigma)\nabla[L]) + 2\left(q_G - \frac{q_{O_2}}{6}\right), \qquad (30) $$

where $D_L$ is the diffusion coefficient (see Table S1 in S1 Document), $q_G$ the glucose consumption rate (Eq 21) and $q_{O_2}$ the oxygen consumption rate (Eq 22). $\sigma$ is the local cell density.

**Extra-cellular matrix.** Components of the ECM are produced by resident cells, and secreted into the ECM via exocytosis [81]. Once secreted, they aggregate with the existing matrix. In the model ECM is produced by cells, diffuses within the interstitial space and degrades with a constant rate

$$ \frac{\partial[\text{ECM}]}{\partial t} = D_{\text{ECM}}\Delta[\text{ECM}] + k_{\text{gen}}^{\text{ECM}}\sigma - k_{\text{deg}}^{\text{ECM}}[\text{ECM}], \qquad (31) $$

where is the density of viable cells, [ECM] is dimensionless and proportional to the ECM concentration divided by a reference concentration which in Col IV stained micrographs corresponds to the highest color intensity. Thus [ECM] ranges from 0 to 1 permitting a direct comparison with the image data.

**Waste material.** Cells dying either by apoptosis or necrosis are finally lysed and release their intra-cellular material into the interstitial space. We assumed some of these molecules inhibiting for cell cycle progression of other cells. The dynamics was modeled by a reaction diffusion equation:

$$ \frac{\partial[W]}{\partial t} = D_W \Delta[W] + k_{gen}^W \sigma_D - k_{upt}^W \sigma[W], \qquad (32) $$

where $\sigma_D$ is the local density of dead cells (no matter if by necrosis or apoptosis) but not yet lysed cells, [W] is the density of waste, $D_W$ is the waste molecule diffusion constant, $k_{gen}^W$ is the rate at which waste molecules leak into the interstitial space, $k_{put}^W$ is the rate at which waste molecules are taken up by viable cells.

## Temporal Evolution & Numerical Treatment

The cellular dynamics was modeled by a master equation describing the time evolution of the probability of the multi-cellular configuration. The molecular dynamics was modeled by a system of deterministic partial differential equations.
**Cell kinetic: Gillespie algorithm.** The master equation is solved numerically by the Stochastic Simulation Algorithm proposed in [82]. The idea is to investigate the expected time interval until the next event (cell migration, cell cycle progression, necrosis, apoptosis, division, ...) and then choose an event randomly, where the probability for each choice reflects its relative weight compared to all possible configuration changes. A detailed description of this part can be found in [83].

**Molecular kinetic: Steady state solution.** As the temporal scale of the cellular dynamics (> hours) and molecular dynamics (< minutes) differ by several orders of magnitude, we could assume that the molecular concentrations reached their steady state very quickly after a change in the cellular configuration. Consequently, we updated the molecular gradients with a constant time step of 0.1h by solving the system of PDEs (Eqs 19, 28, 30, 31 and 32) for their steady state solution. We verified that choosing a smaller time step did not have any significant impact on the model kinetics, and by solving the full time dependent equation for selected parameter combinations, that the time scale at which the concentrations reached their steady state values were short compared to the time scale of cell growth, division and death considered here. To calculate the steady state solutions all the PDEs were solved by a finite difference scheme, with a central difference discretization of the diffusion operator, on the regular lattice used to construct the unstructured lattice populated by the cells (see S1 Document). The coupled equations for glucose and oxygen were solved monolithically, with a semi-implicit iterative scheme. The other PDEs were solved with the same spatial discretization. BiCGSTAB was used to solve all the linear problems.

## Supporting Information

S1 Document. Details about the tumor growth model, model development, additional simulation results and additional information on the image analysis.
(PDF)

## Author Contributions

Conceived and designed the experiments: NJ BM IEVC MM DD. Performed the experiments: BM. Analyzed the data: NJ. Contributed reagents/materials/analysis tools: NJ BM IEVC MM DD. Wrote the paper: NJ BM IEVC DD.

## References

1. Drasdo D, Hohme S. A single-cell-based model of tumor growth in vitro: monolayers and spheroids [Journal Article]. Phys Biol. 2005; 2(3):133-47. doi: 10.1088/1478-3975/2/3/001 PMID: 16224119

2. Schaller G, Meyer-Hermann M. Multicellular tumor spheroid in an off-lattice Voronoi-Delaunay cell model. Physical Review E. 2005; 71(5):051910. doi: 10.1103/PhysRevE.71.051910

3. Pampaloni F, Reynaud EG, Stelzer EH. The third dimension bridges the gap between cell culture and live tissue [Journal Article]. Nat Rev Mol Cell Biol. 2007; 8(10):839-45. doi: 10.1038/nrm2236 PMID: 17684528

4. Breslin S, O'Driscoll L. Three-dimensional cell culture: the missing link in drug discovery [Journal Article]. Drug Discov Today. 2013; 18(5-6):240-9. doi: 10.1016/j.drudis.2012.10.003 PMID: 23073387

5. Tung YC, Hsiao AY, Allen SG, Torisawa YS, Ho M, Takayama S. High-throughput 3D spheroid culture and drug testing using a 384 hanging drop array [Journal Article]. Analyst. 2011; 136(3):473-8. Available from: http://pubs.rsc.org/en/Content/ArticleLanding/2011/AN/c0an00609b. doi: 10.1039/C0AN00609B PMID: 20967331

6. Evensen NA, Li J, Yang J, Yu X, Sampson NS, Zucker S, et al. Development of a high-throughput three-dimensional invasion assay for anti-cancer drug discovery [Journal Article]. PLoS One. 2013; 8(12):e82811. Available from: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3859626/pdf/ pone.0082811.pdf. doi: 10.1371/journal.pone.0082811 PMID: 24349367
7. Friedrich J, Ebner R, Kunz-Schughart LA. Experimental anti-tumor therapy in 3-D: spheroids–old hat or new challenge? [Journal Article]. Int J Radiat Biol. 2007; 83(11–12):849–71. Available from: http://informahealthcare.com/doi/abs/10.1080/09553000701727531. doi: 10.1080/09553000701727531 PMID: 18058370

8. Hirschhaeuser F, Menne H, Dittfeld C, West J, Mueller-Klieser W, Kunz-Schughart LA. Multicellular tumor spheroids: an underestimated tool is catching up again [Journal Article]. J Biotechnol. 2010; 148 (1):3–15. doi: 10.1016/j.jbiotec.2010.01.012 PMID: 20097238

9. Lin RZ, Chang HY. Recent advances in three-dimensional multicellular spheroid culture for biomedical research [Journal Article]. Biotechnol J. 2008; 3(9–10):1172–84. Available from: http://onlinelibrary.wiley.com/doi/10.1002/biot.200700228/abstract. doi: 10.1002/biot.200700228 PMID: 18566957

10. Ernst A, Hofmann S, Ahmadi R, Becker N, Korshunov A, Engel F, et al. Genomic and expression profiling of glioblastoma stem cell-like spheroid cultures identifies novel tumor-relevant genes associated with survival [Journal Article]. Clin Cancer Res. 2009; 15(21):6541–50. Available from: http://clincancerres.aacrjournals.org/content/15/21/6541.full.pdf. doi: 10.1158/1078-0432.CCR-09-0695 PMID: 19861460

11. Kim H, Phung Y, Ho M. Changes in global gene expression associated with 3D structure of tumors: an ex vivo matrix-free mesothelioma spheroid model [Journal Article]. PLoS One. 2012; 7(6):e39556. Available from: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3380922/pdf/ pone.0039556.pdf. doi: 10.1371/journal.pone.0039556 PMID: 22737246

12. Kunz-Schughart LA, Freyer JP, Hofstaedter F, Ebner R. The use of 3-D cultures for high-throughput screening: the multicellular spheroid model [Journal Article]. J Biomol Screen. 2004; 9(4):273–85. Available from: http://jbx.sagepub.com/content/9/4/273.full.pdf. doi: 10.1177/1087057104265040 PMID: 15191644

13. Mehta G, Hsiao AY, Ingram M, Luker GD, Takayama S. Opportunities and challenges for use of tumor spheroids as models to test drug delivery and efficacy [Journal Article]. J Control Release. 2012; 164 (2):192–204. Available from: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3436947/pdf/ nihms378814.pdf. doi: 10.1016/j.jconrel.2012.04.045 PMID: 22613880

14. Freyer JP, Sutherland RM. Regulation of growth saturation and development of necrosis in EMT6/Ro multicellular spheroids by the glucose and oxygen supply. Cancer research. 1986; 46(7):3504–3512. PMID: 3708582

15. Wise SM, Lowengrub JS, Frieboes HB, Cristini V. Three-dimensional multispecies nonlinear tumor growth—I: model and numerical method. Journal of theoretical biology. 2008; 253(3):524–543. doi: 10.1016/j.jtbi.2008.03.027 PMID: 18485374

16. Preziosi L, Tosin A. Multiphase modelling of tumour growth and extracellular matrix interaction: mathematical tools and applications. Journal of mathematical biology. 2009; 58(4–5):625–656. doi: 10.1007/s00285-008-0218-7 PMID: 18853162

17. Basan M, Risler T, Joanny JF, Sastre-Garau X, Prost J. Homeostatic competition drives tumor growth and metastasis nucleation [Journal Article]. HFSP J. 2009; 3(4):265–72. doi: 10.2976/1.3086732 PMID: 20119483

18. Chen Y, Wise SM, Shenoy VB, Lowengrub JS. A stable scheme for a nonlinear, multiphase tumor growth model with an elastic membrane [Journal Article]. Int J Numer Method Biomed Eng. 2014; 30 (7):726–54. doi: 10.1002/cnm.2624 PMID: 24443369

19. Maini P, Deutsch A, Dormann S. Cellular Automaton Modeling of Biological Pattern Formation: Characterization, Applications, and Analysis. Modeling and Simulation in Science, Engineering and Technology. Birkhauser; 2007. Available from: http://books.google.de/books?id=8yRDPFOTDnsC

20. Anderson A, Rejniak K. Single-Cell-Based Models in Biology and Medicine. Mathematics and Biosciences in Interaction. Birkhäuser Basel; 2007. Available from: https://books.google.de/books?id=B9-s6AY1Z-kC. doi: 10.1007/978-3-7643-8123-3

21. Drasdo D, Hoehme S, Block M. On the role of physics in the growth and pattern formation of multi-cellular systems: What can we learn from individual-cell based models? Journal of Statistical Physics. 2007; 128(1–2):287–345. doi: 10.1007/s10955-007-9289-x

22. Jiang Y, Pjesivac-Grbovic J, Cantrell C, Freyer JP. A multiscale model for avascular tumor growth. Biophysical journal. 2005; 89(6):3884–3894. doi: 10.1529/biophysj.105.060640 PMID: 16199495

23. Alarcón T, Byrne HM, Maini PK. A design principle for vascular beds: the effects of complex blood rheology. Microvascular research. 2005; 69(3):156–172. doi: 10.1016/j.mvr.2005.02.002 PMID: 15896358

24. Drasdo D, Kree R, McCaskill J. Monte Carlo approach to tissue-cell populations. Physical Review E. 1995; 52(6):6635. doi: 10.1103/PhysRevE.52.6635

25. Block M, Scholl E, Drasdo D. Classifying the expansion kinetics and critical surface dynamics of growing cell populations [Journal Article]. Phys Rev Lett. 2007; 99(24):248101. doi: 10.1103/PhysRevLett.99.248101 PMID: 18233492
26. Kim Y, Stolarska MA, Othmer HG. A Hybrid Model for Tumor Spheroid Growth in vitro I: Theoretical Development and Early Results [Journal Article]. Mathematical Models and Methods in Applied Sciences. 2007; 17:1773-1798. Available from: http://www.worldscinet.com/m3as/17/17supp01/S0218202507002479.html. doi: 10.1142/S0218202507002479

27. Welter M, Bartha K, Rieger H. Emergent vascular network inhomogeneities and resulting blood flow patterns in a growing tumor [Journal Article]. J Theor Biol. 2008; 250(2):257-80. doi: 10.1016/j.jtbi.2007.09.031 PMID: 17996256

28. Ramis-Conde I, Drasdo D, Anderson AR, Chaplain MA. Modeling the influence of the E-cadherin-beta-catenin pathway in cancer cell invasion: a multiscale approach [Journal Article]. Biophys J. 2008; 95 (1):155-65. Available from: http://www.ncbi.nlm.nih.gov/pubmed/18339758. doi: 10.1529/biophysj.107.114678

29. Ramis-Conde I, Chaplain MA, Anderson AR, Drasdo D. Multi-scale modelling of cancer cell intravasation: the role of cadherins in metastasis [Journal Article]. Phys Biol. 2009; 6(1):016008. Available from: http://www.ncbi.nlm.nih.gov/pubmed/19321920. doi: 10.1088/1478-3975/6/1/016008

30. Szabó A, Merks RM. Cellular potts modeling of tumor growth, tumor invasion, and tumor evolution. Frontiers in oncology. 2013; 3. doi: 10.3389/fonc.2013.00087 PMID: 23596570

31. Drasdo D, Jagiella N, Ramis-Conde I, Vignon-Clementel I, Weens W, et al. Modeling steps from a benign tumor to an invasive cancer: examples of intrinsically multi-scale problems. Cell Mechanics: From Single Scale-Based Models to Multiscale Modeling, Chapman & Hall/CRC. 2010; p. 379-417.

32. Bartha K, Rieger H. Vascular network remodeling via vessel cooption, regression and growth in tumors [Journal Article]. J Theor Biol. 2006; 241(4):903-18. doi: 10.1016/j.jtbi.2006.01.022 PMID: 16545398

33. Lee DS, Rieger H, Bartha K. Flow correlated percolation during vascular remodeling in growing tumors [Journal Article]. Phys Rev Lett. 2006; 96(5):058104. doi: 10.1103/PhysRevLett.96.058104 PMID: 16486998

34. Hoehme S, Drasdo D. Biomechanical and nutrient controls in the growth of mammalian Cell populations. Mathematical Population Studies. 2010; 17(3):166-187. doi: 10.1080/08898480.2010.491032

35. Angelini ED, Clatz O, Mandonnet E, Konukoglu E, Capelle L, Duffau H. Glioma dynamics and computational models: a review of segmentation, registration, and in silico growth algorithms and their clinical applications. Current Medical Imaging Reviews. 2007; 3(4):262-276. doi: 10.2174/157340507782446241

36. Konukoglu E, Clatz O, Menze BH, Stieltjes B, Weber MA, Mandonnet E, et al. Image guided personalization of reaction-diffusion type tumor growth models using modified anisotropic eikonal equations. Medical Imaging, IEEE Transactions on. 2010; 29(1):77-95. doi: 10.1109/TMI.2009.2026413

37. Swanson KR, Bridge C, Murray J, Alvord EC. Virtual and real brain tumors: using mathematical modeling to quantify glioma growth and invasion. Journal of the neurological sciences. 2003; 216(1):1-10. doi: 10.1016/j.jns.2003.06.001 PMID: 14607296

38. Frieboes HB, Smith BR, Chuang YL, Ito K, Roettgers AM, Gambhir SS, et al. An integrated computational/experimental model of lymphoma growth. PLoS computational biology. 2013; 9(3):e1003008. doi: 10.1371/journal.pcbi.1003008 PMID: 23555235

39. Macklin P, Edgerton ME, Cristini V. 10. In: Cristini V, Lowengrub JS, editors. Agent-based cell modeling: application to breast cancer. Cambridge, UK: Cambridge University Press; 2010. p. 206-234. (invited author: P. Macklin). Available from: http://MathCancer.org/Publications.php#macklin10_cup10

40. Macklin P, Edgerton ME, Thompson AM, Cristini V. Patient-calibrated agent-based modelling of ductal carcinoma in situ (DCIS): from microscopic measurements to macroscopic predictions of clinical progression. Journal of theoretical biology. 2012; 301:122-140. doi: 10.1016/j.jtbi.2012.02.002 PMID: 22342935

41. Macklin P, Mumenthaler S, Lowengrub J. Modeling multiscale necrotic and calcified tissue biomechanics in cancer patients: application to ductal carcinoma in situ (DCIS). In: Multiscale Computer Modeling in Biomechanics and Biomedical Engineering. Springer; 2013. p. 349-380.

42. Drasdo D, Hoehme S, Hengstler JG. How predictive quantitative modelling of tissue organisation can inform liver disease pathogenesis. Journal of hepatology. 2014; 61(4):951-956. doi: 10.1016/j.jhep.2014.06.013 PMID: 24950483

43. Drasdo D. Coarse graining in simulated cell populations. Advances in Complex Systems. 2005; 8 (02n03):319-363. Available from: http://www.worldscientific.com/doi/abs/10.1142/S0219525905000440. doi: 10.1142/S0219525905000440

44. Grasl-Kraupp B, Ruttkay-Nedecky B, Koudelka H, Bukowska K, Bursch W, Schulte-Hermann R. In situ detection of fragmented DNA (TUNEL assay) fails to discriminate among apoptosis, necrosis, and autolytic cell death: a cautionary note. Hepatology. 1995; 21(5):1465-1468. PMID: 7737654
45. Lopez-Meraz ML, Niquet J, Wasterlain CG. Distinct caspase pathways mediate necrosis and apoptosis in subpopulations of hippocampal neurons after status epilepticus. Epilepsia. 2010; 51(s3):56-60. doi: 10.1111/j.1528-1167.2010.02611.x PMID: 20618402

46. Higuchi M, Honda T, Proske RJ, Yeh ET. Regulation of reactive oxygen species-induced apoptosis and necrosis by caspase 3-like proteases. Oncogene. 1998; 17(21):2753-2760. doi: 10.1038/sj.onc.1202211 PMID: 9840939

47. Radszuweit M, Block M, Hengstler JG, Scholl E, Drasdo D. Comparing the growth kinetics of cell populations in two and three dimensions [Journal Article]. Phys Rev E Stat Nonlin Soft Matter Phys. 2009; 79(5 Pt 1):051907. doi: 10.1103/PhysRevE.79.051907 PMID: 19518480

48. Drasdo D, Hoehme S. Modeling the impact of granular embedding media, and pulling versus pushing cells on growing cell clones. New Journal of Physics. 2012; 14(5):055025. doi: 10.1088/1367-2630/14/5/055025

49. Kempf H, Bleicher M, Meyer-Hermann M. Spatio-temporal cell dynamics in tumour spheroid irradiation [Journal Article]. The European Physical Journal D. 2010; 60(1):177-193. doi: 10.1140/epjd/e2010-00178-4

50. Galle J, Loeffler M, Drasdo D. Modeling the effect of deregulated proliferation and apoptosis on the growth dynamics of epithelial cell populations in vitro [Journal Article]. Biophys J. 2005; 88(1):62-75. Available from: http://www.ncbi.nlm.nih.gov/pubmed/15475585. doi: 10.1529/biophysj.104.041459 PMID: 15475585

51. Byrne H, Drasdo D. Individual-based and continuum models of growing cell populations: a comparison. Journal of mathematical biology. 2009; 58(4-5):657-687. doi: 10.1007/s00285-008-0212-0 PMID: 18841363

52. Caccavari F, Valdembri D, Sandri C, Bussolino F, Serini G. Integrin signaling and lung cancer. Cell adhesion & migration. 2010; 4(1):124-129. doi: 10.4161/cam.4.1.10976

53. Sanders MA, Basson MD. Collagen IV-dependent ERK activation in human Caco-2 intestinal epithelial cells requires focal adhesion kinase [Journal Article]. J Biol Chem. 2000; 275(48):38040-7. doi: 10.1074/jbc.M003871200 PMID: 10986280

54. Berg J, Tymoczko J, Stryer L. Biochemistry 5th edition. W H Freeman; 2002.

55. Sholl-Franco A, Fragel-Madeira L, Macama AdCC, Linden R, Ventura ALM. ATP controls cell cycle and induces proliferation in the mouse developing retina. International Journal of Developmental Neuroscience. 2010; 28(1):63-73. doi: 10.1016/j.ijdevneu.2009.09.004 PMID: 19799993

56. Vermeulen K, Van Bockstaele DR, Berneman ZN. The cell cycle: a review of regulation, deregulation and therapeutic targets in cancer. Cell proliferation. 2003; 36(3):131-149. doi: 10.1046/j.1365-2184.2003.00266.x PMID: 12814430

57. Mariaule G, Belmont P. Cyclin-Dependent Kinase Inhibitors as Marketed Anticancer Drugs: Where Are We Now? A Short Survey. Molecules. 2014; 19(9):14366-14382. doi: 10.3390/molecules190914366 PMID: 25215591

58. Chen JB, Liu WJ, Che H, Liu J, Sun HY, Li GR. Adenosine-5'-triphosphate up-regulates proliferation of human cardiac fibroblasts. British journal of pharmacology. 2012; 166(3):1140-1150. doi: 10.1111/j.1476-5381.2012.01831.x PMID: 22224416

59. Ozturk SS, Riley MR, Palsson BO. Effects of ammonia and lactate on hybridoma growth, metabolism, and antibody production. Biotechnology and bioengineering. 1992; 39(4):418-431. doi: 10.1002/bit.260390408 PMID: 18600963

60. Lao MS, Toth D. Effects of ammonium and lactate on growth and metabolism of a recombinant Chinese hamster ovary cell culture. Biotechnology progress. 1997; 13(5):688-691. doi: 10.1021/bp9602360 PMID: 9336989

61. Chen C, Byrne H, King J. The influence of growth-induced stress from the surrounding medium on the development of multicell spheroids. Journal of mathematical biology. 2001; 43(3):191-220. doi: 10.1007/s002850100091 PMID: 11681526

62. Helmlinger G, Netti PA, Lichtenbeld HC, Melder RJ, Jain RK. Solid stress inhibits the growth of multicellular tumor spheroids. Nature biotechnology. 1997; 15(8):778-783. doi: 10.1038/nbt0897-778 PMID: 9255794

63. Stolper DA, Revsbech NP, Canfield DE. Aerobic growth at nanomolar oxygen concentrations. Proceedings of the National Academy of Sciences. 2010; 107(44):18755-18760. doi: 10.1073/pnas.1013435107

64. Greenspan H. Models for the growth of a solid tumor by diffusion. Stud Appl Math. 1972; 51(4):317-340. doi: 10.1002/sapm1972514317

65. Hoehme S, Brulport M, Bauer A, Bedawy E, Schormann W, Hermes M, et al. Prediction and validation of cell alignment along microvessels as order principle to restore tissue architecture in liver
regeneration [Journal Article]. Proc Natl Acad Sci U S A. 2010; 107(23):10371-6. doi: 10.1073/pnas.0909374107 PMID: 20484673

66. Andriani F, Garfield J, Fusenig NE, Garlick JA. Basement membrane proteins promote progression of intraepithelial neoplasia in 3-dimensional models of human stratified epithelium. International journal of cancer. 2004; 108(3):348-357. doi: 10.1002/ijc.11525

67. Novaro GM, Aronow HD, Militello MA, Garcia MJ, Sabik EM. Benzocaine-induced methemoglobinemia: experience from a high-volume transesophageal echocardiography laboratory [Journal Article]. J Am Soc Echocardiogr. 2003; 16(2):170-5. doi: 10.1067/mje.2003.5 PMID: 12574744

68. Miao LJ, Huang SF, Sun ZT, Gao ZY, Zhang RX, Liu Y, et al. MiR-449c targets c-Myc and inhibits NSCLC cell progression [Journal Article]. FEBS Lett. 2013; 587(9):1359-65. doi: 10.1016/j.febslet.2013.03.006 PMID: 23507140

69. Warburg E. [Cardiac insufficiency and its treatment] [Journal Article]. Nord Med. 1956; 56(37):1351. WARBURG, E dan Not Available 1956/09/13 Nord Med. 1956 Sep 13;56(37):1351. PMID: 13378710

70. Zu XL, Guppy M. Cancer metabolism: facts, fantasy, and fiction [Journal Article]. Biochemical and Bio-physical Research Communications. 2004; 313(3):459-465. Available from: http://www.sciencedirect.com/science/article/pii/S0006291X0302504X. doi: 10.1016/j.bbrc.2003.11.136 PMID: 14697210

71. Kasinskas RW, Venkatasubramanian R, Forbes NS. Rapid uptake of glucose and lactate, and not hypoxia, induces apoptosis in three-dimensional tumor tissue culture [Journal Article]. Integr Biol (Camb). 2014; 6(4):399-410. doi: 10.1039/c4ib00001c

72. Rodriguez-Enriquez S, Torres-Marquez ME, Moreno-Sanchez R. Substrate oxidation and ATP supply in AS-30D hepatoma cells [Journal Article]. Arch Biochem Biophys. 2000; 375(1):21-30. doi: 10.1006/abbi.1999.1582 PMID: 10683245

73. Olavarria JS, Chico E, Gimenez-Gallego G, Nunez de Castro I. Effect of ammonium ions on the aerobic glycolysis in Ehrlich ascites tumor cells [Journal Article]. Biochimie. 1981; 63(6):469-75. doi: 10.1016/S0300-9084(81)80079-X PMID: 7260172

74. Medina MA, Sanchez-Jimenez F, Marquez FJ, Perez-Rodriguez J, Quesada AR, Nunez de Castro I. Glutamine and glucose as energy substrates for Ehrlich ascites tumour cells [Journal Article]. Biochem Int. 1988; 16(2):339-47. Available from: http://www.ncbi.nlm.nih.gov/pubmed/3365266. PMID: 3365266

75. Guppy M, Kong SE, Niu X, Busfield S, Klinken SP. Method for measuring a comprehensive energy budget in a proliferating cell system over multiple cell cycles [Journal Article]. J Cell Physiol. 1997; 170 (1):1-7. doi: 10.1002/(SICI)1097-4652(199701)170:1%3C1::AID-JCP1%3E3.0.CO;2-S PMID: 9012779

76. Beucher S, Lantuéjoul, C. Use of watersheds in contour detection. In: International workshop on image processing, real-time edge and motion detection; 1979.

77. Vincent L, Soille P. Watersheds in Digital Spaces: An Efficient Algorithm Based on Immersion Simulations [Journal Article]. IEEE Trans Pattern Anal Mach Intell. 1991; 13(6):583-598. doi: 10.1109/34.87344

78. Eden M. A two-dimensional growth process. Dynamics of fractal surfaces. 1961; 4:223-239.

79. Freyer JP, Sutherland RM. A reduction in the in situ rates of oxygen and glucose consumption of cells in EMT6/Ro spheroids during growth [Journal Article]. J Cell Physiol. 1985; 124(3):516-24. doi: 10.1002/jcp.1041240323 PMID: 4044662

80. Casciari JJ, Sotirchos SV, Sutherland RM. Variations in tumor cell growth rates and metabolism with oxygen concentration, glucose concentration, and extracellular pH [Journal Article]. J Cell Physiol. 1992; 151(2):386-94. doi: 10.1002/jcp.1041510220 PMID: 1572910

81. Plopper G. 19. In: Lewin B, Cassimeris L, Lingappa V, G P, editors. The extracellular matrix and cell adhesion. Jones and Bartlett; 2007. p. 821-880.

82. Gillespie DT. Exact stochastic simulation of coupled chemical reactions [Journal Article]. The Journal of Physical Chemistry. 1977; 81(25):2340-2361. doi: 10.1021/j100540a008

83. Alfonso JC, Buttazzo G, Garcia-Archilla B, Herrero MA, Nunez L. Selecting radiotherapy dose distributions by means of constrained optimization problems [Journal Article]. Bull Math Biol. 2014; 76 (5):1017-44. doi: 10.1007/s11538-014-9945-7 PMID: 24599739