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, 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.
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 have now been established as experimental systems for both basic research and high throughput 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), viability promoters (VP), growth inhibitors (GI) and viability inhibitors (VI). In the case of spheroids 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 co-workers [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.
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. 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 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). 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.
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 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−17 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−17 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 contast 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, O2 = [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, [O2] = 0. 28mM, and [G] = 25mM, [O2] = 0. 08mM correctly. 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 (θ) = ∏ i 1 2 π σ i 2 e x p ( (x (θ) i - μ i) 2 2 σ i 2), 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. 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 O2). 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.
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% CO2 in 150cm2 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% CO2 and either normal atmospheric 20% O2 (corresponding to 0. 28mM) or 5% O2 (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. The images acquired (see above) are raw data. They consist of a set of pixels with position (x, y) ∈ N 2 and color intensities for the three color channels red, green and blue defined by I c h a n n e l [ x, y ] ∈ [ 0,1), I c h a n n e l: N 2 → R, (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. The cellular automaton model, which extends on previous work ([43,47]), is defined by a set of rules: 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: ∂ t u = ∇ · (D u (σ) ∇ u) + r u (σ, u), (19) where Du is the molecule diffusion coefficient and ru (σ, u) denotes the reaction term, u ∈ {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 (σ, u) = - q u (u) σ = - q u (u) ∑ k δ (r - r k). (20) 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.
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.
carbohydrate metabolism
cell death
chemical compounds
oxygen
cell cycle and cell division
cell processes
carbohydrates
glucose metabolism
organic compounds
glucose
oxygen metabolism
cellular structures and organelles
extracellular matrix
chemistry
biochemistry
cell biology
organic chemistry
apoptosis
monosaccharides
biology and life sciences
physical sciences
metabolism
chemical elements
Inferring Growth Control Mechanisms in Growing Multi-cellular Spheroids of NSCLC Cells from Spatial-Temporal Image Data