Blood Vessel Tortuosity Selects against Evolution of Agressive Tumor Cells in Confined Tissue Environments: a Modeling Approach
Andr\'as Szab\'o, Roeland M. H. Merks
TL;DR
This study uses a computational model to explore how blood vessel tortuosity and blood flow obstruction influence tumor evolution, revealing that vascular instabilities can suppress aggressive tumor cells and affect tumor development.
Contribution
The paper introduces a novel cell-based model linking blood flow obstruction to tumor cell evolution, highlighting the impact of vascular instabilities on tumor progression and phenotype selection.
Findings
Blood flow obstruction can lead to tumor cell extinction.
Vascular instabilities influence tumor development stages.
Aggressive phenotypes may be suppressed by blood vessel tortuosity.
Abstract
Cancer is a disease of cellular regulation, often initiated by genetic mutation within cells, and leading to a heterogeneous cell population within tissues. In the competition for nutrients and growth space within the tumors the phenotype of each cell determines its success. Selection in this process is imposed by both the microenvironment (neighboring cells, extracellular matrix, and diffusing substances), and the whole of the organism through for example the blood supply. In this view, the development of tumor cells is in close interaction with their increasingly changing environment: the more cells can change, the more their environment will change. Furthermore, instabilities are also introduced on the organism level: blood supply can be blocked by increased tissue pressure or the tortuosity of the tumor-neovascular vessels. This coupling between cell, microenvironment, and organism…
| Parameter | Notation | Initial | Step | Minimum | Maximum |
| name | value | size | value | value | |
| Incompressibilty | 25 | 0.2 | 0 | 20 | |
| Division volume | 50 | 1 | 10 | 200 | |
| Cell adhesion | 99 | 0.1 | 0 | 100 | |
| Matrix adhesion | 49.5 | 0.1 | 0 | 50 | |
| Glucose chemotaxis | 0 | 0.1 | -20 | 20 | |
| Oxygen chemotaxis | 0 | 0.1 | -20 | 20 | |
| Lactate chemotaxis | 0 | 0.1 | -20 | 20 | |
| Growth (see Eq. 13) | 250 | 25 | 0 | 2500 | |
| HIF threshold | 0.1 | 10-3 | 10-4 | 10 | |
| ROS threshold | 0.5 | 0.01 | 10-4 | 100 |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Blood Vessel Tortuosity Selects against Evolution of Agressive Tumor Cells in Confined Tissue Environments: a Modeling Approach
András Szabó1*¤ and Roeland M. H. Merks1,2
1 Life Sciences Group, Centrum Wiskunde & Informatica, Amsterdam, The Netherlands 2 Mathematical Institute, Leiden University, Leiden, The Netherlands
¤Current Address: Department of Cell and Developmental Biology, UCL, London, United Kingdom
Abstract
Cancer is a disease of cellular regulation, often initiated by genetic mutation within cells, and leading to a heterogeneous cell population within tissues. In the competition for nutrients and growth space within the tumors the phenotype of each cell determines its success. Selection in this process is imposed by both the microenvironment (neighboring cells, extracellular matrix, and diffusing substances), and the whole of the organism through for example the blood supply. In this view, the development of tumor cells is in close interaction with their increasingly changing environment: the more cells can change, the more their environment will change. Furthermore, instabilities are also introduced on the organism level: blood supply can be blocked by increased tissue pressure or the tortuosity of the tumor-neovascular vessels. This coupling between cell, microenvironment, and organism results in behavior that is hard to predict. Here we introduce a cell-based computational model to study the effect of blood flow obstruction on the micro-evolution of cells within a cancerous tissue. We demonstrate that stages of tumor development emerge naturally, without the need for sequential mutation of specific genes. Secondly, we show that instabilities in blood supply can impact the overall development of tumors and lead to the extinction of the dominant aggressive phenotype, showing a clear distinction between the fitness at the cell level and survival of the population. This provides new insights into potential side effects of recent tumor vasculature renormalization approaches.
Author Summary
Multicellular organisms control their cells to facilitate higher level function of the whole organism. In tumors, this control is lost and cells are allowed to enhance their fitness by, for example, increased proliferation. Tumor cells continue to change their behavior through accumulating mutations, leading to a complex and highly heterogeneous structure. Several computational studies have investigated the emergent structures of such mutating group of cells and led to the recognition that cellular heterogeneity within tumors is essential to explain the observed morphologies. Most of these studies have considered a limited number of possible cell phenotypes, an isolated tumor cell population, unlimited growth space around the tumor, often with an unexhaustable source of nutrients. Here we introduce a modelling approach that takes into account the limited growth space around the tumor, localized nutrient sources, cellular metabolism, and mutation in a continuous phenotype space. The model reproduces the Warburg effect due to the limited nutrient supply leading to an irreversible switch in cellular metabolism, and, consistently with previous models, exhibits stages of development together with a natural selection for a rapidly growing phenotype. This phenotype locally emerges in stable environments, but when nutrient supply becomes erratic, these show less resilience and are outcompeted by slow growers.
Introduction
Cancer is a disease of multicellular regulation, in which malfunctioning cells can break free of homeostatic regulations imposed by the host environment [1]. One of the main characteristics of cancer is the increased proliferation and mutation of cancerous cells due to malfunctioning control of growth and proliferation [1]. As these behavioral changes typically originate from mutations in the cells’ genetic material, excessively proliferating cells accumulate further alterations, leading to a possible amplification of malignancies. Traditional studies of altered cell traits primarily focus on genetic mutations, but neglect the multicellular nature and genetic variety of tumors.
Tumor heterogeneity has been demonstrated experimentally and is an active field of research [2, 3, 4]. Neutral mutations may accumulate and contribute to intratumor heterogeneity [5]. An intermediate level of heterogeneity is correlated with low survival probability [6]. Heterogeneity may even promote the collapse of tumor development by inducing a clone population that supports and enhances the growth of other clones in mice [7]. Marusyk and colleagues [7] claim that these supporter clones may be outcompeted by the more aggressive subpopulation, leading to the disappearance of the supporters and to the consequent collapse of the tumor. Heterogeneity questions the validity of previous whole–tumor analyses, as “the most abundant cell type might not necessarily predict the properties of mixed populations” [8], and emphasizes the need for more detailed approach.
Initial phases of tumor development are increasingly thought to give rise to a Darwinian process [1, 9], where individual cells compete for growth space and nutrients. Selection is imposed by the microenvironment, a highly complex entity spanning several ranges in size from the endocrine regulation of the whole body down to the extracellular matrix (ECM) and neighboring cells. During cancer development this environment is changed due to changes in the cellular component, and due to the tissue and organism level reactions to the tumor.
Intrinsic coupling between neighboring cells forms the basis of the plasticity-reciprocity model [10]: as cells alter their behavior (plasticity), their contribution to the local environment changes through for example ECM remodelling, nutrient uptake, or adhesion molecule expression. In turn, the changed environment imposes an altered selection on the cells (reciprocity), creating a feedback between cells and their local microenvironment. This cascading change in behavior is reminiscent of the behavioral changes associated with stages of cancer development [11] that was later described by the accumulation of mutations through which cells become increasingly malignant [12]. Although a strict sequence of mutations was not found, a general pattern was observed in the majority of cases [12], linking the macroscopic stages to the cell-level changes.
Computational modeling is an excellent tool for exploring, studying, and understanding such complex systems, because it provides complete control over assumptions and reveal the consequent behavior emerging from them. This allows the dissection of complex interactions and exploration of experimentally challenging cases. A variety of models have been applied to study cancer.
Using a population level description, Basanta and co-workers studied how the combination of tumor treatments, p53 cancer vaccine and chemotherapy, can be optimized to yield the best results [13]. In a similar study, Sreemati Datta and coworkers [14] model tumor development with the inclusion of evolving mutation rates and show that the balance between inducing driver mutations and mutation rates plays a key role in tumor growth: at high mutation rates, genetic instability may counter tumor progression. Combined with close experimental verification, Marusyk and colleagues [7] used a similar model to suggest that interactions among clones may lead to an overall collapse of tumor development. This approach is able to incorporate evolutionary games to help cope with the development of treatment resistance. Such space-free models assume that all cells within the tumor are able to interact with all other cells directly, and are unable to deal with intra-tumor spatial heterogeneities.
Studies incorporating this spatial aspect have mostly worked with cellular automata (CA) models. For example, Gerlee and colleagues were able to explain the ‘go or grow’ hypothesis through the emergence of haptotaxis in their CA model [15]. In a recent study, Waclaw and colleagues [16] used a CA model to show that cell motility together with cell turnover may prevent intratumor heterogeneity. Of particular interest is the study of Anderson and colleagues [17], where the evolution of a growing cell population and the effect of a heterogeneous environment is explored. They represent the tumor environment as a distribution of ECM molecules that together with oxygen serve as nutrient after degradation. Cell evolution is modeled in the phenotype by selecting a set of cellular parameters (matrix degradation rate, proliferation rate, etc.) from a number of predefined phenotypes upon cell division. The new phenotype was either selected randomly or according to a predefined sequence progressing towards more aggressive behavior. The authors found that cells evolved into a similar, aggressive phenotype when applying the random mutation scheme. Heterogeneity in the population was found to give rise to irregular tumor surface, whereas environmental heterogeneity (heterogeneity in the ECM distribution) reduced population heterogeneity and favored the most aggressive cell types. Low concentration of oxygen was also found to reduce population heterogeneity and promote invasive finger formation; applying two bursts of oxygen in these simulations led to the segregation of the mixed population. Similar models have been used to show emergent progression of phenotypes related to hypoxia, glycolysis, and acid-resistance, by including neuronal networks in cells [18], or angiogenesis by introducing blood vessels in hypoxic regions [19]. Enderling and co-workers [20] introduced the cancer stem cell (CSC) hypothesis in a similar model by incorporating cells with unlimited proliferation potential (CSCs) and cells with limited proliferation potential. They found that cell death induced by tumor therapy could lead to a more aggressive, proliferative tumor, as the CSCs were no longer competing for space after the treatment. Using a combination model of phenotype evolution and the CSC hypothesis, Sottoriva and colleagues showed that the presence of CSCs in the model tumors led to more invasive tumor morphology [21]. The CA model further allowed exploring efficient drug delivery through the neovasculature [22, 23], the role of spatial arrangement of the vasculature in radiation therapy [24] and even the effect of the different 2D representations of the 3D vasculature [25].
Using a more detailed model allows to explore the evolution of many other aspects of cell behavior, such as cell flexibility, cell adhesion, or cell shape. One such model used in tumor modeling is the cellular Potts model [26, 27]. It has been used to describe, for example, the effect of nutrient limitation on tumor growth morphology [28], or to compare the emergence of distinct developmental stages in terms of morphology and growth in vascular versus avascular tumors [29]. Using a heterogeneous evolutionary model of CSCs, a recent study reported the emergence of spatial stratification of tumors based on the evolution of adhesion molecule expression [30]. A further line of models explore the effect of oxygenation of tumors [31] in the light of chemotherapy [32]. An alternative way of introducing more detail is a spatial continuum representation of the cells. One example of this is the phase field method which has been used in combination with discrete modelling elements to investigate tumor growth morphology and its interaction with the vasculature [33, 34].
Most of the above models applies a sharp distinction between the tumor and the healthy tissue. While the tumor is described in a cell-based detail, the environment is usually presented as a continuum. These models implicitly assume that stromal cells surrounding the tumor do not participate in the development of the tumor, and that these cells are fundamentally different from the cancerous tissue. The transition between the stromal and cancerous cell states is neglected. A recent exception from this is the work of Powathil and colleagues [35] where the effect of irradiation on a small set of “healthy” bystander cells is investigated. As a result of irradiation and induced signals, stromal cells may apoptose or may be converted to the tumor cell type in the model. However, in this and almost all of the above models the phenotypes available for evolving cells are restricted to a small, discrete subset of all possibilities. A continuous range of cell states could reflect a more biological picture including more subtle, for example epigenomic, changes. Furthermore, the nutrients are typically modeled as a single component, however, one signature of cancerous cells is the reduced oxygen consumption and increased glucose uptake. For this a more detailed nutrient description is required. Another less explored point of interest is the change in the temporal behavior of the larger environment, the nutrient supply. Initial tumors experience a stable blood supply in healthy tissues, that presents them with a relatively constant environment. In contrast, angiogenic tumors at a later stage develop neovasculature that is tortuous and leaky, presenting a fluctuating, unstable environment for the cells [36, 37]. How does this temporal variation affect the population behavior? Would the same aggressive phenotype dominate the population as in the stable environment, or would it fit less than the less aggressive and stable phenotype?
Here we present a cellular Potts model of a closely packed, mutating cell population representing an epithelial tissue within an organism (Fig. 1). Mutation in the model allows cellular behaviors to vary continuously in a wide range of phenotype space, therefore evolution is governed by selection emerging naturally from the limitations of space, glucose, and oxygen. Cellular metabolism is modeled as a mixture of aerobic glycolysis and respiration. Small scale changes in the microenvironment are represented by the changes of the immediate cellular neighborhood, both by mutations and cell rearrangements. We model large scale environmental fluctuations through the fluctuating activity of spatially fixed nutrient sources, which represent the cross-section of blood vessels.
Constructed in this way our model includes both the plasticity-reciprocity model of Friedl and Alexander [10], as well as the large scale fluctuations imposed by the host. In the following sections we analyse the behavior of the model with and without mutating cells. We show that our model reproduces the Warburg shift, and exhibits stages of development similar to those observed in previous studies (such as [17, 18, 19]). In our model cells undergo clonal expansion, hypoxia, followed by starvation, with the development of segregated populations around blood vessels. The spatial differentiation of cell populations is somewhat similar to the spatial diversity in real tumors as described by Alfarouk et al. [38]. Whereas Alfarouk and colleagues describe two main habitat zones concentrically surrounding the blood vessel, we observe only one of the zones with high proliferation rates and a robust cellular outflow from near the nutrient source. Finally, our results indicate that the dominant aggressive phenotype is more sensitive to fluctuations in the environment than the ones maintaining a stable phenotype without mutation.
Results
Cellular Potts model of a homeostatic tissue
To investigate the above questions, we model a monolayer of cells using a modified cellular Potts model (CPM) based on the CompuCell3D implementation [39] which can be obtained from http://www.compucell3D.org. Customized code for the simulations and example parameter and initial condition files can be found in S1 File.. In the following we give an overview of the model; for more detail see the Methods section.
Cells in the CPM are represented as confluent domains on a lattice on which an integer at every position indicates which cell is occupying the location ; cell-free areas are designated by . Cell movement results from a series of elementary steps in which an attempt is made to copy at a randomly selected location to one of its randomly selected neighboring location . This attempt is accepted with a probability based on a Hamiltonian goal function that defines cell dynamics (Eqs. 1, 2). is usually defined such that cells maintain a controlled size, perform amoeboid-like cell movement, and may exhibit adhesion or contact-repulsion. A time step in the model is defined as the Monte Carlo Step (MCS) consisting of elementary steps where is the total number of lattice sites in the model. In our model we apply the usual calibration by relating 1 MCS to 1 minute real time, and 1 lattice site to 2 m. Diffusion of soluble substances are simulated on a lattice identical to the cellular-lattice. This calibration relates the simulated tissue area to , and diffusion coefficients of simulated nutrients to realistic values (glucose and lactate: ; oxygen: ) [40].
We implemented a metabolism whereby cells consume glucose and oxygen from their environment and use a mixture of lactic acid fermentation and cellular respiration. Cells metabolize oxygen and glucose to generate an abstract cellular energy that is used for their maintenance and growth. The amount of energy required for cells is controlled by the expression levels of glucose transporters in the cell membrane which is determined by an intracellular growth signal parameter ( for cell at time ). The mode of metabolism is determined by an internal hypoxia inducible factor ( for cell at time ), that is controlled through the oxygen levels inside and around the cell and the amount of intracellular reactive oxygen species (ROS)(Eq. 15). This factor determines the ratio of respiratory and fermentative modes of cellular metabolism. In our model the speed of energy production, and hence cell growth, is independent of the mode of metabolism in order to avoid a selection bias towards the faster metabolic mode.
If a cell reaches a pre-defined doubling size it divides, hence cell cycle time is determined by cell metabolism. Cell death may occur either due to starvation or age. If a cell generates less energy from metabolism than is required for maintenance, it converts the necessary amount of its cell mass into energy (catabolism). Once the cell mass is exhausted, the cell is considered dead and is taken out of the simulation. Cells are also killed in the simulation with a 0.1% probability after each MCS to maintain cell turnover.
Glucose and oxygen are supplied by a separate set of designated immobilized cells that play the role of blood capillaries (Fig. 1a). To allow the temporal control of nutrient supply, capillaries can be in an active or a blocked state. Nutrient levels are kept at a fixed concentration in the blood stream, that is inside the capillaries, when the capillaries are active. When a capillary is blocked, nutrient levels are kept at zero within the capillaries. The activity of the vessels are changed with a probability (nutrient switching probability) after each MCS. A high switching probability leads to a stable supply while a low switching probability results in an inconsistent supply, mimicking blocked or tortuous vessels without affecting the average activity time of the vessels. Lactate produced by cells as a waste from lactic acid fermentation is cleared out of the system by active capillary cells.
For more detail on the model see the Methods section.
We start by first verifying that the model is capable of simulating a sustainable homeostatic tisse. Therefore we studied the behavior of a healthy tissue in the absence of mutation and stable nutrient supplies. Fig. 2a shows the model setup. We consider a two-dimensional square lattice, corresponding with a slice of tissue of 400 m 400 m (40 40 cells, 200 200 lattice sites), containing four blood vessels arranged in a square formation (Fig.1a). To achieve constant nutrient supply, the probability of a vessel to be blocked or unblocked in every MCS is 0.5. In this case each vessel switches between active (depicted in white) and blocked (depicted in gray) states rapidly. As the nutrients diffuse away from the source, this rapid switching results in a continuous supply of nutrients (Fig. 2c). The color of tissue cells in Fig. 2a indicates the intracellular pressure, defined as the difference between target volume (biomass) and actual volume of the cell (). This measure differs from the extracellular pressure as it includes any contribution of the contractile actin cortex surrounding biological cells. Note that the pressure within the population is distributed without any specific pattern.
Fig. 2b shows the dynamics of the tissue for MCS corresponding to approximately 70 days. The number of cells fluctuates around a constant value throughout the simulation (Fig. 2b), and is sufficient to cover the whole system: the ratio of the cell-covered region over the total simulation area is close to one (Fig. 2b). As the average intracellular pressure does not increase over the course of the simulation (Fig. 2b), cell growth and proliferation are kept in balance with the basal metabolism and the constant cell turnover. Excessive growth is prevented by the lack of growth space through a negative feedback between the intracellular pressure and cell growth (Eq. 18). Nutrient levels remain constant during the simulations and cells remain respiratory as indicated by the absence of lactate (Fig. 2d).
Model of tissue micro-evolution
Next we investigated how a cancerous tissue would behave in our model. Cancerous tissues are characterized by large number of mutations, chromosomal rearrangements and changes in gene expression level, all of which may result in phenotypic changes of the cells [41, 42, 1, 43]. To mimic such changes, we allowed the set of 10 assigned phenotypic properties of cancerous cells to change upon division (Fig. 1b), including the division volume or adhesion parameters (see Table 1 and Methods). After cell division, the daughter cells inherit the phenotypes of their parents with some small mutations. Every parameter is allowed to change with a fixed probability (mutation rate for parameter ) and independently of one another. The change in parameter is drawn from a normally distributed random variable with a standard deviation of , that is: . The parameters are allowed to change freely within a pre-defined range (Table1) with reflective boundary conditions. This allows an unbiased parameter to uniformly explore the available range (Fig. 1b inset).
To determine how a single mutating cell would perturb the homeostasis of the above tissue, we inserted a cell with a mutation potential either near (Fig 3a, red cell) or far from (Fig 3b) the nutrient source. The single mutating cell persisted in 28 and 26 simulations out of 100 for the two different initiation positions. When the cell persisted, it expanded within the first 5000 MCSs and eventually colonised the population completely (Fig 3a, b).
In order to study the internal dynamics of the tumor, we will only focus on the case where the mutating cell persists and has colonized the population; therefore we will initiate our simulations with populations where all cells are allowed to mutate. Fig. 3c-f shows the behavior of the model with mutating cell populations, with 10% mutation rate. In comparison with the non-mutating populations, the intracellular pressure is higher and exhibits a wider variation across neighbors (Fig. 3c). This shows that cells overcome the initial growth control mechanism. The number of cells initially increases, reaches a peak, and then declines to approximately half of the peak value, well below the initial numbers, where the population size stabilizes (Fig. 3d). These changes in cell numbers are followed by the intracellular pressure as well: in the expansion phase the pressure increases, but before the peak in cell number it sharply declines. After the population size settles to a lower value, the pressure settles to an approximately constant value. The full coverage of the tissue drops to approximately half coverage, showing that the cancerous cells cannot maintain a complete monolayer. Nutrients are depleted further from the blood vessels, resulting in a shortage of oxygen in the distant regions (Fig. 3e), shortly followed by the depletion of glucose (Fig. 3f). The depletion of oxygen triggers the cells to switch to fermentation, resulting in an increase in lactate. This switch accelerates the depletion of glucose, causing the decline in population size. After the population size is reduced, oxygen levels return to the same level as in the non-mutating populations, while cells still rely on fermentation as can be seen from the maintained lactate levels (Fig. 3f).
Emergent stages of development
Previous studies have reported distinct stages of development including hypoxia, glycolysis, or acid-resistance [17, 18, 19]. However, in these studies the evolution occurred in isolation from the stromal tissues and vasculature either using a limited set of phenotypic behaviors potentially constraining the degree of freedom of the evolutionary trajectories or with immobile cells. Therefore, we first asked if stages of tumor progression also occurred in our less constrained model. Fig. 4 shows the behavior of our model with the nutrient concentrations and cell numbers averaged from 10 independent simulations with stable vasculature (switching probability = 0.5) and 10% mutation rate. Based on this, we identified distinct stages in our model: expansion (1), hypoxia (2), starvation (3). Insets show cell configurations characteristic of the three stages, color scale on the insets indicates the oxygen concentrations. These stages emerge as a result of an interplay between the cells and their environment, as in the proposed plasticity-reciprocity hypothesis of Friedl and Alexander [10]. Stages in our model relate to: 1. Conditioning the environment; 2. A reaction to the environmental change in the behavior of cells (new phenotypes emerging, old phenotypes disappearing); 3. New environment created by the new population. Remarkably, this shows that despite the much larger number and freedom of mutating parameters in our model, we still find the same phenomena of emergent stages.
In the first stage of our model the population expands by cell growth and division. A high intracellular growth signal is selected for in the population, favoring fast growing cells (Fig. 4c). This parameter evolves much faster than any other of the 10 mutating parameters (see also Table 1) while most of the other parameters do not exhibit such a strong and clear drift by the end of the simulations (Fig. 4d). Indeed, the overall behavior of the model did not change qualitatively in simulations where only or and the chemotactic parameters are allowed to evolve (Supporting Information). Since nutrients in the environment are not limiting due to the assumed prior homeostatis, these cells simply outgrow the slower ones, creating patches of high growth (Fig. 4e). This leads to the expansion of fast growers, and as a result, cells from newer generations appear in clumps of growth hot-spots (Fig. 4f). At this stage expansion can occur at any position in the population since nutrients are available at any location.
In the second stage of our model the population turns hypoxic. As the number of cells grows, the intracellular pressure increases rapidly in the tightly packed tissue until about t = 5000 MCS (see Fig. 3d). At this time oxygen is depleted at areas further away from the source (see Fig. 4 middle inset at t=6000 MCS). Fast growing cells (high , Fig. 4g) are unable to fuel their increased metabolic need through oxidative respiration and turn hypoxic (Fig. 4h). These cells further increase their glucose uptake (Fig. 4i) due to the HIF1- GLUT signaling pathway in our model (Eq. 15), and start the production of lactate.
Finally, glucose is gradually depleted at regions far from the sources as a result of the elevated glucose consumption rate of fast growing cells. In the depleted areas cells die out, and with them the cell population is gradually decreased (Fig. 3b,d). The only cells remaining are around the vessels, that eventually hijack the source (Fig. 4j, k). These cells continue to compete as in stage 1, since the change in the environment near the vessels is minimal, and keep increasing their internal growth signal from generation to generation (Fig. 4j). Increasing intracellular pressure near vessels (Fig. 4k) exerted by neighboring cells and counteracts growth.
Changes in nutrient levels indicate that our model selects for cells exhibiting the well-known Warburg effect, whereby cells metabolize glucose through glycolysis even in the presence of oxygen (aerobic glycolysis) [44]. Cells in our model initially shift to aerobic glycolysis to support their metabolic need escalated through competition, shown by the increasing levels on intracellular hypoxia and ROS (Fig. 4l). This results in an increase in extracellular oxygen (Fig.4a). Despite the availbility of oxygen, cells are unable to revert to a more efficient full respiratory metabolism due to production of ROS which stabilizes HIF1 and limits the amount of metabolic flux through respiration (Eq. 15), thus keeping it in a state of hypoxia in our model (Fig. 4l). Nevertheless, cells do consume oxygen but it is significantly lower than glucose uptake (Fig. 4m).
Taken together, these results show that our model exhibits different stages of development similar to previously published studies. Remarkably, this progression emerges in spite of an almost completely unrestricted evolution of a large number of phenotypic parameters. Tumors in this model are initialized at random positions, but due to the explicit representation of localized nutrient sources, we show that they occupy the vicinity of blood vessels at later stages. This is enhanced by the more realistic representation of cells in the CPM where cell shape and compressibility allow cell rearrangements within the packed tissue as opposed to the more rigid CA models exploring progression [17, 18, 19]. Secondly, we show that our model selects for cells exhibiting the Warburg effect despite the lack of growth advantage of fermenting cells.
Higher mutation rate speeds up transition between stages
To test if the stages of tumor progression depend on the phenotypic mutation rates, we simulated the model for a series of mutation rates. Whereas the non-mutating population keeps a constant size, all mutating populations exhibit an initial increase in cell numbers (Fig. 5a). In highly mutating populations (5% and 10%) this increase is followed by a drop in cell numbers. This drop is observed later in populations with 5% mutation rate, and population decrease is just starting at the end of the simulations in populations with 1% mutation rate. Note that the repetitions reproduce the behavior fairly well, suggesting the robustness of the system. Therefore we suggest that similar stages occur at lower mutation rates, and the time needed for reaching each stage depends on the mutation rate. This is supported by the changes in nutrient levels in the simulations, which react faster to change than the total number of cells (Fig. 5b-f). In healthy, non-mutating populations, nutrient levels and population size is stabilized (Fig. 2b). In mutating populations, cells evolve a higher metabolic demand in parallel with increased proliferation. This results in an increase in population size and decrease in oxygen levels, followed by a decrease in glucose levels. As oxygen becomes sparse, cells turn hypoxic and switch from oxydative phosphorylation to aerobic glycolysis, resulting in an increase in lactate levels. We have found the same behavior in populations with different mutational probabilities ranging from 0.1% up to 10% (Fig. 5b-f), or higher (Supporting Information a-h). Again, a lower mutational probability only delayed the changes in the nutrient levels. We did not find a qualitative difference in the progression in our simulations at different mutation rates, showing that the emergent order of stages is robust in this system.
To explore the structure of the population at different mutation rates in the system, we analysed the distribution of cells in the space of normalized mutating parameters. After subtraction of the initial parameter values and normalizing with mutational step-size, the ten-dimensional parameter space of the population was reduced to the three most prominently changing axes within each population using principal component analysis (see Methods). Fig. 5g shows one example population at the final time point of the simulations for mutation rates 1%, 10%, and 20% with each dot representing a cell. At low mutation rate (1%) the population splits up into well-defined clones which are more spread and less well-defined at higher mutation rates. The first three principal axes contain most of the information about the shape of the population, as can be seen by the normalized weights (eigenvalues) of these components (Fig. 5h). The composition of principal axes at the final time point of the simulations shown on Fig. 5g differ (Fig. 5i) with the growth signal , cell rigidity , and glucose chemotaxis playing an important role at 1%. At 10% mutation rate the lactate chemotaxis parameter plays a distinct role in segregating the population phenotypes, while at the segregation is less obvious (Fig. 5g) and is driven mainly by parameters , doubling volume , and adhesions (). While the composition of the main axes varies across different simulation repeats with the same mutation rate, the populations are nevertheless well characterized by the first three components in all cases (Supporting Information i, j).
To better understand population structure in the simulations, we categorized the cells at each time point in the 10-D phenotype space using hierarchical clustering (Methods). We measured the displacement of each cluster as the Euclidean distance between the point of origin and its center of mass and its spread as the mean distance of points of the cluster from the cluster’s center of mass. As expected, we observed that populations with higher mutation rates reach further from the origin by the end of simulations; these clusters are more spread and less dense than clusters in populations with lower mutation rates (Fig. 5j, k). Considering the whole time course of the simulation, the cluster analysis reveals that as the population explores the phenotype space, clusters from the highly mutating populations first tend to spread out and dilute more than the clusters in populations from lower mutation rates (Fig. 5l, m). These results show that the population starts to dilute much faster in phenotype space at high mutation rates, but without affecting the progression of stages apparent from the nutrient levels.
Faster growing cell population is less robust in the face of external nutrient fluctuations
Next, we tested how feedback from a larger spatial organizational level, through the nutrient supply in our case, would affect populations of different mutation rates. In a healthy tissue, nutrient supply is relatively constant. The main source of fluctuations are the slow daily change according to the circadian rhythm, and the relatively fast blood pulse. In cancerous tissues the vasculature is remodeled through tumor vasculogenesis, resulting in tortuous and leaky vessels [36, 37]. As these vessels are less reliable, here we assume that they dysfunction from time to time, for example by becoming temporarily blocked.
To model vessel tortuosity, we introduced a blocking probability for the vessels. An open blood vessel will be blocked with a probability at every time step, and a blocked vessel will be opened with the same probability. A high blocking probability () corresponds to a healthy situation, where the resulting fast switching of the vessel is smoothened out by nutrient diffusion. A lower blocking probability introduces longer periods of nutrient deprivation but also longer periods of nutrient supply. On average these systems receive the same amount of nutrients, but in different dosage.
Healthy, non-mutating populations in our simulations survived over a wide range of nutrient fluctuations. The average number of cells from 10 simulation repeats showed that these populations keep roughly the same size at different blocking probabilities, as shown on Fig. 6a. At the extreme blocking probability , only 2 out of 10 populations died out.
In comparison, mutating cell populations are unable to tolerate blocking probabilities lower than , irrespective of their mutation rate (Fig. 6a). Note that at high , populations with low mutation rates have an increased population size at the end of the simulations ( MCS), compared to the healthy population. This increase results from the initial stage of progression, as these populations only reach the first stage (expansion) by the end of the simulations. However, this advantage disappears as decreases.
The observed reduction in population size due to decreasing could work in two ways: either by killing cells through starvation, or by speeding up the progression of stages. In the previous section we showed that a higher mutation rate speeds up the progression of the population and thus results in a reduced population size (Fig. 5). This might eventually lead to extinction. If the populations under fluctuating nutrient supply go through the same stages of progression as the ones in the stable environment, the environmental indicators used in Fig. 5 (levels of glucose, oxygen, lactate) are unsuitable, as these might change in simulations with different . Instead, we focus on the intracellular evolution of traits, that are not altered directly in these experiments. The first trait to be selected for is the intracellular growth signal of the cells () that exhibits a run-away dynamics (Fig. 4c). If the blocking probability accelerates the progression through the stages, it should increase the selection pressure on the intracellular growth signal as well. Contrary to this expectation, we found that the trend in the average value of the intracellular growth signal remains approximately the same in simulations across different values, and even slightly decreases at lower (Fig. 6b). Similarly, other cellular measures (such as hypoxia, ROS, or pressure), or cellular parameters (such as the chemotaxis parameters) show the same behaviors irrespective of (Supporting Information). Therefore, longer nutrient fluctuations do not accelerate the evolution of the population, and the reduction in population size is not a result of the acceleration of the same evolutionary dynamics. Instead we conclude that as cells deplete the nutrients in the environment due to their increased consumption, the chance for survival in systems with longer fluctuations is reduced.
Next we altered during the time of simulation runs. We tested how the population reacts if the blood vessels become increasingly tortuous, starting from a healthy state (fast switching) progressing to a tortuous vasculature (slow switching). Blocking probability in these simulations is decreased gradually in the simulations, following a geometric progression with an initial value of and ratio of . Once the progression reaches at MCS the blocking probability is not decreased further (). Similar to the stable system, the mutating populations (mutation rate ) are initially driven into the high consumer state by cell-cell competition (Fig. 6c solid red line showing growth signal parameter ). When the fluctuation probability reaches the order of ( MCS), the populations start to die out (Fig. 6c red dashed line showing number of surviving populations), similar to the case of static low blocking probabilities. Note however, that non-mutating populaions (Fig. 6c blue) are able to survive increasing fluctuations in the nutrient supply. This shows that a changing nutrient supply does not necessarily influence the direct competition among cells.
Inconsistency of nutrients may emerge from the disfunctional tissues of the emergent tumor occluding vasculature. To represent this feedback, we examined how the population behaves when the consistency of nutrient supply is related to the amount of cellular coverage in the tissue. We created a feedback between the density of the tissue (measured as tissue surface coverage , with ) and the fluctuating source. For a fully populated tissue we kept the nutrient switching probability high, , providing a smooth nutrient supply, and decreased it linearly with the cell density to model the variability in nutrient supply. Thus: . In these simulations the mutating population persists at approximately the same level as in the healthy case (Fig. 6d). The emergent tissue coverage yields an approximate fluctuation probability of , sufficiently high to support mutating populations. While the relationship between nutrient supply stability and living cell density is experimentally unclear, our results suggest that a linear relationship between stability and density is insufficient to drive the agressive cells to extincition, similar to what is expected of an expanding pathological tumor.
The effect of low nutrient levels with two consecutive surges of high oxygen levels has been shown to affect selection in a model of evolutionary tumor growth [17]. Here we used localized nutrient sources and stochastic supply, rather than two deterministic and uniform pulses. We show that in our system, inconsistent nutrient supply does not increase the speed of cellular evolution, however, it does reduce the viability of the unstable cell populations. We found similar results when inconsistency is considered progressively increasing or proporitional to the tissue coverage. Our results suggest that, while consistent nutrient supply promotes cell-level selection of fast growing cells, inconsistent nutrient supplies that put a larger demand on the tissue exert selection at the tissue scale and provide higher chances of survival for populations with cellular quiescence.
Clonal expansion and heterogeneity as a consequence of differential growth
In a simulated healthy tissue devoid of mutants, cell growth is independent of spatial localization, growth is limited intrinsically by a constant growth signal balancing spatial confinement. After this intrinsic limitation is lifted through uncontrolled mutations and competition, growth becomes clustered in growth hot-spots (Fig. 4e). These spots are populated by overly proliferative cells, resulting in clonal expansion: Fig. 7a shows configurations at two time points in the same simulation color-coded for descendants. At a later stage, when cells deplete nutrients, proximity to the sources dictates the growth rate (Fig. 4j). In the resulting environment cells closer to the source grow faster and divide, while cells further away starve and die. This differential growth gives rise to a directed cell movement from the vicinity of the sources to the depleted areas, apparent from the short cell trajectories shown in Fig. 7b. To demonstrate that the segregation patterns were indeed linked to the nutrient sources and are not the result of finite system size, we performed simulations with randomly scattered blood vessels instead of the regular distribution, similar to previous studies [24, 31]. The example trajectory plot on Fig. 7c shows that the outflow of cells is correlated with the nutrient positions.
The population becomes segregated, as the probability of a cell invading a neighboring region is diminished. As a result of this differential growth cells near the source take over the vicinity of the vessel and spread their phenotype in this region (Fig. 7d). While cells within each segregated part of the population have highly similar phenotype, the different parts evolve independently of one another. This is shown by the distribution of phenotypic traits over time in each population (Fig. 7e). This gives rise to independently evolving quasi-species-like populations within the tissue, and might be analogous to the observed heterogeneity in tumors.
In simulations with low but constant vessel blocking probabilities cells around a blocked source are able to leave their region and migrate to another active source, as exemplified by the trajectories in Fig. 7e. The ability to detect and to move to a neighboring active source is crucial for the cells for survival, and therefore is expected to be selected for. In our model, cells can achieve this by evolving chemotaxis towards glucose or oxygen, or chemotaxis away from lactate. Indeed, chemotaxis parameters are the second most affected parameters throughout the evolution of the population (Fig.4c, d). Interestingly, chemorepulsion by lactate is strongly selected for although lactate in our model does not have any direct metabolic effect on the cells. This acquired property helps to orient the cells towards the nutrient sources better than oxygen, which becomes ubiquitous with more shallow gradients, or glucose, which does not diffuse far from the nutrient sources due to elevated uptake. The sum effect of these motility parameters can be expressed with a combined chemotaxis parameter defined as: . Fig. 7g and h show the evolution of this combined chemotaxis parameter in populations from different vessel blocking probabilities and mutation rate of 10%. Indeed in all of the conditions the combined chemotaxis is selected for. In simulations with lower blocking probabilities the combined chemotaxis has increasingly higher fluctuations due to the random selection introduced by the blocking and opening of the vessels, often leading to extinction as well.
Discussion
In this study we introduced a computational model to explore the evolution of tissue cells under specific conditions. Selection in the model acts through physical constraints: limited growth space, and limited nutrient diffusion. Selected cellular phenotype traits, such as cell rigidity or growth signal, are allowed to change (“mutate”) upon cell division in any direction with a probability within a given range.
Progression of phenotypes
We show that a directional evolution emerges from random movement in phenotype space as the result of cell competition, driving cells from a healthy (homeostatic) state to a more aggressively expansive phenotype. This is consistent with previous findings of Anderson and colleagues [17], who showed that a similar drift towards aggressive phenotypes emerges if cells are allowed to mutate randomly into one of a 100 predefined discreet phenotypes. In contrast to these abrupt changes in behavior, our model only allows small changes in mutation. This choice lends more persistence to the clones in the population, since more mutation events are required to diverge from an existing clone. Faster growing cells are selected in our model, which then go on to colonize the population by means of clonal expansion.
Progression of phenotypes has been observed previously in other models of tumor evolution, where the authors also considered the toxic effect of acidification due to glycolysis [18, 19]. In [18], cells first develop the ability to survive in hypoxic environments by lowering their apoptitic response threshold to low oxygen levels. This is followed by a metabolic switch to glycolysis and finally the emergence of the acid-resistant phenotype. In contrast, acid-resistance was proposed to emerge first followed by evolution of glycolysis in another model where nutrient sources were represented as point sources [19]. The distinct stages exhibited in our model (clonal expansion, oxygen depletion, and starvation) are accompanied by a sudden increase in glycolytic activity which is then moderated onto a stable medium level (Fig. 4l). In our model the population turns hypoxic (Fig. 4l) before developing a distinct chemo-repulsive response (Fig. 4c). We found the progression robust against changing the mutation rate, however, at higher mutation rates the population was able to explore the phenotype space faster and in more dispersed clusters (Fig. 5).
Although lactate in our model does not affect cells in our model, cells gradually develop a chemorepulsion away from lactate (a negative on Fig. 4c, d). Lactate accumulates at regions where high consumers deplete glucose, therefore cells use lactate chemorepulsion as a compass to navigate away from the high-consumer niche towards a more supportive environment. This novel feature in our model highlights how phenotypes that show no apparent advantage at the cell level could be selected for based on the altered micro-environmental conditions.
Minor differences in proliferation lead to clonal expansion
Local expansion in tumor growth models has recently been described by Waclaw and colleagues [16]. Their study focused on heterogeneity in passanger mutations while the probabilistically occurring driver mutations set the proliferation advantage. Cell motility there acts to blur intratumor heterogeneity by allowing proliferating cells to invade ‘cell-free’ areas where further proliferation is not inhibited by other cancer cells. While this model is able to reproduce the clonal expansion features of tumors, it neglects the potential inhibitory effect of spatial constraint produced by the surrounding stromal tissue. Proliferation diversity in tumors was shown to be essential to explain experimentally observed tumor morphology by an earlier study using the cancer stem cell hypothesis [20] which interprets the tumor as a conglomerate of self-metastases. The emergent clonal expansion patterns in our model (Fig. 7a) are reminiscent of these self-metastases and therefore suggest that slight differences of proliferation rates within the population are sufficient to generate these patterns, as opposed to the sharp distinction between cancer stem cells and non-stem cells.
Spatial confinement and localized nutrient sources
Previous models of tumor evolution typically neglect the confining effect of the surrounding healthy tissue although spatial confinement can play an important role when studying the effects of treatment recovery in a tumor with a cancer stem cell population [21]. When a large portion of cells is killed by therapy, the internal cells have access to growth space and are able to regrow the tumor. In other words: space limitation keeps the growth inhibited. To mimic this situation here we restrict the growing tumor to a confined space, making them a model of an in vivo system. We note that simulations in our system do not develop a necrotic core as in other models, as dead cells simply shrink and disappear from the simulations without inducing any signaling response from neighboring cells or without increasing spatial confinement. Therefore our system focuses on live tumor cell evolution where necrotic cells may be considered as extruded from the simulated monolayer region into the underlying necrotic core, removed by the immune response and/or drained by the lymphatics.
The source of nutrients in previous models is typically considered uniform, or is provided in the environment in a dispersed fashion. The source of nutrients are the blood vessels in our model, similar to the more recent studies of Shirinifard and colleagues [29], Powathil and coworkers [31, 32, 35], or using a hybrid cellular automaton model [45]. Two recent studies showed that arrangement of vessels and their 2D representation affects tissue oxygenation and may alter the outcome of a simulated radiotherapy [25, 24]. Nutrient diffusion and utilization create a spatial gradient in the supplies which is translated into a differential growth pattern in our system. Differential growth creates a collective flow of cells outwards from the sources (Fig. 7b, c, f). This outward flow is counter-acted by the emergent chemotactic tendency of cells to move closer to the source (Fig. 4c, d). Chemotaxis in a changing environment provides additional advantage in locating active nutrient sources (Fig. 7f).
Blood vessels in our model are immobilized entities, vascular remodeling is not included. Inconsistency in nutrient supply is implemented as stochastic switching of blood vessel activity, but without feedback from the pressure or hypoxia in the tissue as in the phase field model of Yan and colleagues [34]. Simulations where the blocking probability is a function of time (Fig. 6c) or tissue coverage (Fig. 6d) showed that the progression of stages is not altered in our system. A higher number of connected loops have been shown to emerge in 3D models of tumor angiogenesis than in 2D, predicting that complete blockage of circulation is very rare [23, 22]. Here we studied local vessel blockage which would still have a local effect on our system, even if at a reduced frequency.
Inconsistent nutrient supply
Previously Anderson and colleagues showed that surges of oxygen induce diversification into a population of simulated cells evolving under low oxygen levels [17]. A repeated, second surge was shown to induce phenotypic segregation in these simulations [17]. Here we studied how inconsistency of the nutrient supply affects the population using stochastic (rather than deterministic) switching. In contrast to Anderson’s study, cells do not receive a constant low oxygen supply in our model and are therefore prone to extinction. We found that populations with increasing mutation rate are increasingly sensitive to nutrient fluctuations (Fig. 6a), although the evolutionary trend at the cell level remained largely unaffected (Figs. 6b, Supporting Information). In our model the aggressive phenotype analogous to cancer is a natural consequence of selection on the cell-level and random mutation. When selection pressure is applied on the tissue level (in the form of fluctuating nutrient supplies) the overall fitness (or survival) of the cancerous population proves to be lower than that of the healthy population, potentially due to depletion of ambient nutrient resources as a form of competition. Cancerous cells create an insecure environment that is more sensitive to stress coming from outside the cell population.
This theoretical exploratory study opens questions in how to best approach renormalization of the tumor vasculature. In the model, healthy cells are able to withstand the irregularities in nutrient supply better than cancerous cells. Based on this observation, it is tempting to speculate that irregularities of the tortuous tumor vasculature might serve a similar role in real tumor development. It is important to bear in mind that these results are based on a rudimentary model of tumor development. In addition to the simplifications discussed above, cell-cycle regualtion is overly simplified in our model as opposed to the study of Powathil and colleagues [31] where it is in the focus of the study; cells in our model lack an explicit way to store surplus energy as opposed to the cumulative health-factor of Swat and co-workers [30] representing the cells’ tolerance against starvation. Due to its 2D nature, our model is unable to account for the out of plane transport of nutrients or cells. However, using our quasi-2D model allows the study of several processes taking place in epithelial tissues, where most of the tumors arise. Importantly, the implemented cellular metabolism is overly simplified to enable the exploration of the system. After the foundation of this model framework is set out, it could be expanded with more detailed intracellular metabolic networks, for example, using spatial dynamic flux-balance analysis [46]. These multiscale models will lead to a tighter integration of computational and experimental work similar to the “symbiotic” approach described in a recent angiogenic sprouting study of Cruys and colleagues [47].
Further future work includes the more thorough exploration of different mutation rates for different traits. One of the ultimate goals of computational systems biology is constructing a virtual tissue in order to predict efficient treatment of various diseases. This can be achieved by adding homeostatic mechanisms to the tissues such as contact inhibition of growth, or introducing an internal energy storage as in the study of Swat and colleagues [30]. Inclusion of a dynamic angiogenesis model (e.g. [22, 48, 49, 29, 50]) might involve exploring the roles of different blood vessel placements, or implementing angiogenesis models already available in the same platform. Finally, the system enables the measurement of fitness at different spatial levels and in different (micro- and macro-) environments which could serve as a basis for exploring evolutionary trade-offs in using computational simulations.
Conclusion
In this study we introduce an unbiased evolutionary approach to studying the evolution of interacting tissue cells. The model includes localized source of nutrients (oxygen, glucose) and sink for intermediate metabolites (represented by lactate in our model), and a simplified cellular metabolism including glycolysis and respiration. Cells in the model are spatially confined and no explicit distinction is made between stromal or tumor cells. The model exhibits distinct stages of development with an emergent evolutionary drift towards rapid growth, high glucose consumption, and hypoxia. This is accompanied by a Warburg-effect, whereby cells become unable to return to respiratory metabolism even at high oxygen levels. The simulated tumor exhibits clonal expansion and eventually gives rise to similar phenotypes around each nutrient source which then evolve independently later on. Finally, we found that the emergent rapid growing population is highly sensitive to intermittent nutrient supply, such as caused by leaky tortous blood vessels.
Methods
Details of computational model
Cell volume, adhesion and chemotaxis
The CPM was introduced in the early 1990’s [26, 51] and has evolved into a widely used computational model, applied to various phenomena in developmental biology, cell biology, and tumor biology. Cells in the model are represented as domains on a regular square lattice with an integer spin, , on each lattice site that uniquely identifies the biological cell occupying that site, with identifying the extracellular region (medium). Cell movement is described as a series of elementary steps, in which an attempt is made to copy of a randomly selected site into an adjacent site . The attempt is accepted with probability 1 if it decreases the value of a globally defined Hamiltonian function, , or with Boltzmann probability if it would increase the value of :
[TABLE]
Here is the change in the Hamiltonian due to the attempted copy, and parametrizes the noise introduced by the active motion of cells. Immobilized cells are implemented by rejecting any copy from or into the space occupied by them (). The Hamiltonian function consists of the sum of functions that describe distinct biological features:
[TABLE]
The first term is a volume constraint responsible for maintaining a controlled cell volume which we implemented in a novel way as:
[TABLE]
Here is the volume and is the target volume of cell and describes cell compressibility. Our choice of the volume term differs from the originally proposed volume term [26] as we divide with the target volume of the cell. Using the original formulation results in a synchronization of cell divisions in a compressed monolayer of cells. Using our novel formulation, the daughter cells after a cell division are less compressible than the mother, and thus a cell division does not trigger a wave of divisions by allowing the neighbors to expand more. Note that in a homogeneous population of cells (i.e., , with a constant) this formulation is equivalent to the original one.
The second term represents chemotaxis and is used to calculate a chemotactic bias for the change in the Hamiltonian as introduced by Savill and Hogeweg [52]:
[TABLE]
The function represents the concentration of oxygen () glucose () or lactate () at lattice site , and is the corresponding chemotactic coefficient. Diffusing chemicals are modeled with partial-differential equations, which are solved numerically on a grid matching with the grid of the CPM (see below).
For modeling cell-cell adhesion, the interactions are separated into surface-tension related adhesion and adhesion-molecule related adhesion [53], and are described by the last term of Eq. 2 as:
[TABLE]
where
[TABLE]
and if the sites , belong to the same cell, and otherwise. Above the simplified notation and is used for brevity. is the surface tension related adhesion coefficient between cell-types and present at positions and . As is typically positive, cells tend to minimize their surface area with other cells or the medium, making the adhesion term equivalent to surface tension [51]. Three cell types are used in the simulations: (stromal) cells, endothelial cells, and (cell-free) medium; these will be annotated as: , respectively. The second term of in Eq. 6 describes the adhesion molecule related part of adhesion. The matrix defines the adhesion coefficient between adhesion molecule type and , and represents the density of molecule on cell . Therefore, the maximum amount of adhesion related to adhesion molecules and between cells and is determined by . In our model, we chose the following specific adhesions between the existing cell types:
[TABLE]
Here CAM, MAM and MAS represent general “cell adhesion molecules”, “matrix adhesion molecules” and “matrix adhesion sites”. CAM and MAM are expressed in cells but not the medium. MAS are expressed on medium (matrix) on which cells are allowed to adhere. In our simulations all the values of listed above have the value , the rest of the matrix is set to [math].
The usual time measure in the CPM is the Monte Carlo Step (MCS). One MCS is defined as copy attempts, where is the number of lattice sites in the grid. This choice ensures that on average all sites are updated in every MCS, decoupling the system size and the number of copy-attempts needed to update the whole configuration.
Metabolism
We implemented the following cellular metabolism system into the model.
Cells use glucose and oxygen to produce internal energy that is used for the production of biomass and for supporting basal metabolism. In this section we describe how a cell produces energy from nutrients.
Metabolites in the model are represented by the concentration fields of glucose, oxygen, and lactate, which diffuse and decay as:
[TABLE]
Here is the diffusion coefficient and is the decay constant of the respective species and describes the uptake or secretion of the chemical species according to boundary conditions (in ECs) and cellular metabolism. To ensure the desired trafficking of metabolites, we implemented cellular uptake of glucose and oxygen (sink) and secretion of lactate (source) while cells metabolize as described below. To simulate the diffusion process, the above fields are represented on the same lattice used for representing cells with periodic boundary conditions.
Fields in functioning (active) blood vessel cells (ECs) are set as: , to represent the supply of glucose and oxygen and the removal of lactate. Non-functional (inactive) ECs do not alter nutrient levels. We implemented nutrient switching as follows: at every MCS each EC switches between active and inactive states with a switching probability . Setting results in a frequent switching between active and inactive states providing a reliable steady nutrient supply. Lowering increases the length of inactive and active states, leading to a more unreliable nutrient supply. This way can alter the nutrient supply stability. A stable nutrient supply represents a normal vasculature, whereas tortuous neo-vascular vessels are blocked (inactivated) and unblocked (activated) with a higher variation. Here we relate the amount of tortuosity to the frequency of fluctuations in blood vessel function. Activation (unblocking) and inactivation (blocking) of ECs occur with the same probability, therefore simulations with different blocking frequencies have the same average open and blocking times of blood vessels, and hence the same total amount of nutrients over time.
We implemented two different modes in which cells are able to use the nutrients, mimicking aerobic glycolysis and respiration of biological cells. Briefly, glucose is first degraded into two pyruvates which generates two ATPs. The pyruvate might be converted enzymatically to lactate (glycolysis) without further energy gain, or the cell might degrade it using oxydative phosphorylation (respiration) to gain number of ATP molecules per glucose molecule ().
In the model, cells take up glucose through GLUT transporters using a Michaelis-Menten-like dynamics with cell- and time-specific constants: for cell at time the maximum rate is denoted as and the Michaelis-Menten constant is denoted as . The uptake of glucose from the environment at position per unit time is given by:
[TABLE]
Again, . Here is the extracellular glucose concentration at point and time , is the number of glucose transporters (in arbitrary units) on the surface of cell , and is the cell’s surface (assuming homogeneous GLUT distribution). We implemented the total amount of glucose taken up by cell in each MCS as the sum of uptakes over the cell surface:
[TABLE]
where is the total amount of glucose at the surface of cell ().
A part of the glucose, denoted by for cell at time , is used only in glycolysis to produce 2 lactate molecules and 2 ATP molecules per glucose molecule (). The rest of the glucose () is directed to oxydative phosphorylation after glycolysis to produce ATP molecules per 1 glucose molecule in total. Respiration requires 6 oxygen molecules per 1 glucose molecule, and produces reactive oxygen species (ROS, in cell at time ) as a side-product. Therefore in our implementation the uptake of oxygen is determined by the uptake of glucose and the ratio :
[TABLE]
at every point of the surface of cell at time . Here we use the experimental results of Subczynski and colleagues [54], who showed that oxygen permeability cannot represent a constraint to cellular respiration. Thus, we assume that the internal and external oxygen concentrations equilibrate much faster than the metabolic processes that use internal oxygen. We implemented the energy gain of cell at time as the sum of energy produced by fermentation and respiration, minus the cost of basal metabolism:
[TABLE]
where represents the basal metabolic cost for the cell during a time interval , and is considered to be a fixed constant in the model. Therefore, the variable describes the state of hypoxia in the cell, similar to the transcription factor HIF1-. If , the cell is using all the glucose in the oxidative phosphorylation cycle. If , respiration is completely blocked and the cell changes into a state of complete fermentation in which all glucose is turned over into lactate. In order to compensate for the reduced efficiency in energy production in a controlled manner, when the cell is fermenting, , we assumed the number of glucose transporters is upregulated as well (consistent with experimental observations [55]):
[TABLE]
Here is the factor in the GLUT concentration that is independent of the mode of cell metabolism, and sets the intrinsic growth signal of cell at time . This choice of function describes a controlled compensation for the hypoxia, since the energy gained by cell in a MCS is independent of the mode of metabolism, . Substituting from Eqs. 12, 10, and 13:
[TABLE]
The amount of is calculated based on , the average available glucose concentration over the whole cell surface (before consumption), as the amount of is controlled by the cell as a whole. We set a maximum surface density limit on the GLUT molecules, by limiting parameter to a fixed value.
The hypoxia inducible factor, HIF1-, is induced if the level of oxygen is low in the cell: under normal circumstances PHDs (prolyl hydroxylase) combine with oxygen and inactivate HIF1- [56]. When oxygen is low, the degradation of HIF1- is slower, allowing an increase in HIF1- levels, leading to the reduction of respiration and therefore oxygen use. Furthermore, HIF1- is stabilised by the reactive oxygen species (ROS) produced in the mitochondrial respiratory cycle in biological cells, creating a negative feedback to downregulate respiration if too much ROS is present [56]. To model these dependencies, we implement as a function of intracellular ROS () and average oxygen level in the cell using a combination of Hill functions:
[TABLE]
Here and are parameters describing the threshold and sensitivity of the switch to the levels of ROS () and average oxygen level () in cell at time . Note that the value of in the cell is calculated from values of ROS and oxygen within the cell at the previous time step . This time delay is natural, as the signaling is not expected to instantly alter the biochemical state of the cell. Finally, the level of ROS is increased by the amount of glucose processed in mitochondrial respiration in the previous time step, and the previous values of ROS in the cell:
[TABLE]
where is the decay constant of the ROS inside the cell.
Cell growth
In every time step, the cell is allowed to use its energy to produce biomass, by increasing its target volume as:
[TABLE]
with:
[TABLE]
where is a growth rate parameter, describing the efficiency of anabolism: biomass produced per unit of energy per unit of time, and is the energy lost for the basal metabolism per time step (). The second term expresses a tendency of a model cell to reduce its intracellular density fluctuations: If a cell is compressed, its growth is decreased. The third term expresses the effect of ROS () on cell growth. High levels of ROS are toxic for the cell (interferes with protein folding), forcing the model cell to reduce in size. The function is the Heaviside step-function and is a threshold for ROS levels that a cell can tolerate without damage. If cells are not allowed to grow, for example due to contact inhibition of growth, then the target volume is not increased, but is allowed to decrease.
Irrespective of whether a cell grows or not, the energy of the cell is depleted leaving the cell with no internal energy reserve. If the produced energy in a time step is smaller than the metabolic need of the cell (), the cell’s target volume will be reduced by and its energy is set to zero. This way we include catabolism in our model cells, which allows them to use their target volume as a type of internal energy storage.
Cells divide if their size () is around division size . The probability of allowing cell division is defined by a sigmoid function:
[TABLE]
Here is a parameter setting the width of the fuzzy threshold, and is set to the 5% of the division size: . Upon division, the mother cell divides, and its biomass () is shared equally between the two daughter cells.
To maintain cell turnover, cells are selected for apoptosis at a constant probability (0.1% per cell per MCS). When a cell is selected for apoptosis, its target volume is decreased with a constant rate until it reaches zero. At that point the cell is removed from the simulation and its volume is converted into extracellular space.
Evolution in the model
We implemented evolution in the model as follows. To mimic plasticity, a set of cellular parameters are tested independently for mutation in daughter cells. Such a parameter is allowed to change with a probability (mutation rate ). The new value for the parameter is drawn from a normal distribution centered at the previous value, and with a standard deviation determined by the characteristic step size parameter of the mutation for each parameter: , where is drawn from the normal distribution . The parameter values are bound in within a pre-set parameter range, with reflective boundaries (see Fig. 1b). Constructed this way, the cells perform a random walk in the parameter space through mutations as they progress from generation to generation.
The cellular parameters that are allowed to undergo such mutation are controlling the physical properties of the cells (compressibility through , division size ), adhesion properties (), chemotaxis sensitivities (), growth (through ), and metabolism (sensitivity thresholds ). Values for these parameters are summarized in Table 1.
Initial conditions
Cells are initialized in a monolayer with an initial volume and target volume of 25 lattice sites on a lattice of and four endothelial cells (Fig. 1a). In the initial regime of the simulations, nutrients are allowed to diffuse into the system to obtain a natural distribution resulting in high glucose and oxygen and low lactate levels (Fig. 1b). During this equilibrating time cells are allowed to metabolize, but cannot grow, shrink, move, or mutate. This state represents a stable, homeostatic tissue with sufficient nutrient availability. When the temporal changes in diffusing nutrients is less than 5% in a time interval of 100 MCS, the initial regime is closed, cells are released and the simulation is started. Nutrient fields at the beginning of the initial regime (Fig. 1b) are initialized using pre-generated concentration distributions to expedite equilibration. These initial concentration fields are generated by simulating a cell population in the initial regime starting with zero concentrations: , , . The concentration fields are saved when the total concentration levels remain within 5% over a 100 MCS iteration period: for . Note that no lactate is present as all cells are respiratory in the initial regime. All mutating parameters in all cells are initialized with the same value, after which all cells undergo a mutation attempt to provide an initial heterogeneity to the population.
Adhesion parameters are set to , and for the three region (“cell”) types: stromal cells (c), endothelial cells (EC), and cell-free areas (m). To allow the full range of interactions for cells, we fix the adhesion molecule density values of the extracellular region and the endothelial cells as:
[TABLE]
With this choice the effective values are allowed to change between [math] and for , and . For interaction with the medium, is allowed to change between [math] and .
Diffusion parameters for glucose, oxygen, and lactate were set to and following the approximation of Jiang and coworkers [40], and decay is neglected for all of these chemical species (). Decay rate of ROS is set to a constant per MCS. The amount of ATP produced from 1 glucose molecule through respiration was chosen as . This approximation is a theoretical upper limit for the process, in reality this number is expected to be lower. However, due to the compensation with the number of glucose transporters (see Eq. 13) the exact value of this parameter is not expected to change the behavior of the system.
Parameter analysis
To analyse the population behavior over the simulations we analyzed the evolving parameters in the following way. For every parameter of cell at time we calculated a normalized parameter as where is the characteristic step size of the parameter. This way the normalized parameters reflect their distance from their origin in terms of mutational step-size. The distribution of cell phenotypes in this space was analyzed at time by finding the principal components of the 10-dimensional set of points for all and using singular value decomposition from scientific python (SciPy). The axes corresponding to the first three largest eigenvalues were selected as the main components of the cloud.
Clustering in the normalized phenotype space was performed using hierarchical clustering of SciPy with the Ward method and Euclidean distances. To distinguish clusters we established a cutoff cophenetic distance of 100 manually by evaluating a set of selected dendograms and distribution of clusters plotted on the first three main axes. Displacement of the clusters is calculated as the (10D) Euclidean distance of the center of mass of the cluster from its initial point of origin. Spread of clusters was calculated as the mean distance of points in the cluster from its center of mass. Density of clusters was calculated by dividing the spread by the number of points in the cluster.
Supporting Information
S1 File.
Script and code for running simulations of the model using the CompuCell3D framework (ZIP). Simulations were run using a CompuCell3D installation available from the Indiana University Bloomington and the Biocomplexity Institute (www.compucell3d.org) [39]. Parameters for the simulations are set in XML scripts (simulation.xml). Simulations require initial configuration files for cells and diffusible concentration fields (PIFtemplate.pif, FieldGLUtemplate.txt, FieldOXYtemplate.txt) specified in the XML file, and a set of customized plugins (called ’steppables’ and ’plugins’) to be installed in the CompuCell3D ‘DeveloperZone’. Simulation outputs are a set of files for each sampling time point that contain cell-based information (values of mutating parameters and other cell-level indicators), concentration field data, and cell configuration data.
Acknowledgments
The authors thank Indiana University Bloomington and the Biocomplexity Institute for providing the CompuCell3D modeling environment (www.compucell3d.org) [39]. This work was carried out on the Dutch national e-infrastructure with the support of SURF Cooperative (www.surfsara.nl).
This work was cofinanced by the Netherlands Consortium for Systems Biology (NCSB), which is part of the Netherlands Genomics Initiative/Netherlands Organization for Scientific Research (NWO). The work is also part of the research programme “Innovational Research Incentives Scheme Vidi Cross-divisional 2010 ALW” with project number 864.10.009, which is (partly) financed by the Netherlands Organisation for Scientific Research (NWO). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 11. Hanahan D, Weinberg RA. Hallmarks of cancer: the next generation. Cell. 2011;144(5):646–74. doi:10.1016/j.cell.2011.02.013.
- 22. Gerlinger M, Rowan AJ, Horswell S, Larkin J, Endesfelder D, Gronroos E, et al. Intratumor heterogeneity and branched evolution revealed by multiregion sequencing. The New England journal of medicine. 2012;366(10):883–92. doi:10.1056/NEJ Moa 1113205.
- 33. Sottoriva A, Spiteri I, Piccirillo SGM, Touloumis A, Collins VP, Marioni JC, et al. Intratumor heterogeneity in human glioblastoma reflects cancer evolutionary dynamics. Proceedings of the National Academy of Sciences of the United States of America. 2013;110(10):4009–14. doi:10.1073/pnas.1219747110.
- 44. Yuan Y, Failmezger H, Rueda OM, Ali HR, Gräf S, Chin SF, et al. Quantitative image analysis of cellular heterogeneity in breast tumors complements genomic profiling. Science translational medicine. 2012;4(157):157ra 143. doi:10.1126/scitranslmed.3004330.
- 55. Williams MJ, Werner B, Barnes CP, Graham TA, Sottoriva A. Identification of neutral tumor evolution across cancer types. Nature Genetics. 2016;48(August 2015). doi:10.1038/ng.3489.
- 66. Andor N, Graham TA, Jansen M, Xia LC, Aktipis CA, Petritsch C, et al. Pan-cancer analysis of the extent and consequences of intratumor heterogeneity. Nature Medicine. 2015;22(1):105–113. doi:10.1038/nm.3984.
- 77. Marusyk A, Tabassum DP, Altrock PM, Almendro V, Michor F, Polyak K. Non-cell-autonomous driving of tumour growth supports sub-clonal heterogeneity. Nature. 2014;514(7520):54–58. doi:10.1038/nature 13556.
- 88. Marusyk A, Almendro V, Polyak K. Intra-tumour heterogeneity: a looking glass for cancer? Nature reviews Cancer. 2012;12(5):323–34. doi:10.1038/nrc 3261.