Microenvironmental Variables Must Influence Intrinsic Phenotypic Parameters of Cancer Stem Cells to Affect Tumourigenicity

Since the discovery of tumour initiating cells (TICs) in solid tumours, studies focussing on their role in cancer initiation and progression have abounded. The biological interrogation of these cells continues to yield volumes of information on their pro-tumourigenic behaviour, but actionable generalised conclusions have been scarce. Further, new information suggesting a dependence of tumour composition and growth on the microenvironment has yet to be studied theoretically. To address this point, we created a hybrid, discrete/continuous computational cellular automaton model of a generalised stem-cell driven tissue with a simple microenvironment. Using the model we explored the phenotypic traits inherent to the tumour initiating cells and the effect of the microenvironment on tissue growth. We identify the regions in phenotype parameter space where TICs are able to cause a disruption in homeostasis, leading to tissue overgrowth and tumour maintenance. As our parameters and model are non-specific, they could apply to any tissue TIC and do not assume specific genetic mutations. Targeting these phenotypic traits could represent a generalizable therapeutic strategy across cancer types. Further, we find that the microenvironmental variable does not strongly affect the outcomes, suggesting a need for direct feedback from the microenvironment onto stem-cell behaviour in future modelling endeavours.


Introduction
Heterogeneity among cancer cells within the same patient contributes to tumour growth and evolution. A subpopulation of tumour cells, called Tumour Initiating cells (TICs), or cancer stem cells, has recently been shown to be highly tumourigenic in xenograft models and have some properties of normal stem cells. Evidence continues to emerge that TICs can drive tumour growth and recurrence in many cancers, including, but not limited to, brain [1], breast [2] and colon [3]. These tumour types can be broadly classed as hierarchical as they have been posited to have a hierarchical organisation similar but not identical to nonneoplastic stem-cell (SC) driven tissues. In these hierarchical tumors, TICs can differentiate to produce non-TIC cancer cells or self-renew to promote tumor maintenance. As TICs have been demonstrated to be resistant to a wide variety of therapies including radiation and chemotherapy, the TIC hypothesis has important implications for patient treatments [4]. Specifically, the effect of current strategies on the tumor cell hierarchy should be defined, and TIC specific therapies are likely to provide strong benefit for cancer patients.
In a simplified view of the tumour cell hierarchy, TICs can divide symmetrically or asymmetrically to produce two TIC daughters or a TIC daughter and a more differentiated progeny [5,6], respectively. More differentiated TIC progeny which still have the capability of cell division and are similar to transient amplifying cells (TACs) in the standard stem-cell model and are capable of several rounds of their own symmetric division before the amplified population then differentiates into terminally differentiated cells (TDs) which are incapable of further division. This mode of division and differentiation, which we will call the Hierarchical Model (HM) is schematized in Figure 1.
In the HM, there are a number of cellular behaviours that govern the system. In this study, we choose to study three: the rate of symmetric versus asymmetric division of the stem cells (a), the number of rounds of amplification that transient amplifying cell can undergo before terminal differentiation (b), and the relative lifespan of a terminally differentiated cell (c). While it is a simplification of reality to study only these three parameters and leave out others (for example: differing proliferation rates for the different cell types [7] or the differing metabolic demands of stem vs. non-stem cells [8]) rigorous quantification of these parameters has been extremely difficult to pin down experimentally and so the majority of the work to describe them has been in silico. Most germane to the loss of homeostasis is the work by Enderling et al. [9] which showed the changes to the size of a mutated tissue (tumour) as they varied the number of rounds of amplification of TACs. Other recent work attempting to quantify the rate of symmetric to asymmetric division in putative glioma stem cells was presented by Lathia et al. [10], who showed that this rate can change depending on the presence or absence of growth factors, suggesting yet another method by which a tissue can lose or maintain homeostasis: in reaction to microenvironmental change. A critical limitation of in vivo lineage tracing performed to date is an inability to determine the impact of microenvironmental heterogeneity on TIC symmetric division.
While the HM appears to be quite straight forward, there is growing evidence of complexity to be further incorporated into the model. There are likely to be differences in the extent of TIC maintenance or the ability of tumour cells to move toward a TIC state. TICs appear to reside in distinct niches suggesting there may be differences in the biology of these cells, but defining differences in TICs is limited by cell isolation and tumour initiation methods. Prospective isolation of TICs relies on surface markers, including CD133, CD151 and CD24 which can be transient in nature [11], due to modulation by the tumour microenvironment [12] or methods of isolation [13]. Characterisation of these sorted cells then requires functional assays including in vitro and in vivo limiting dilution assays [14]. In this formulation, each stem can undergo two types of division, either symmetric (with probability a) or asymmetric (with probability 1{a). Each subsequently generated transient amplifying cell (TAC) can then undergo a certain number (b) of round of amplification before differentiating into a terminally differentiated cell (TD) which will live for a certain amount of time before dying (c timesteps). It is these three parameters, which we assume are intrinsic to a given stem cell, which we explore in this paper. doi:10.1371/journal.pcbi.1003433.g001

Author Summary
In this paper, we present a mathematical/computational model of a tumour growing according to the canonical cancer stem-cell hypothesis with a simplified microenvironment. We explore the parameters of this model and find good agreement between our model and other theoretical models in terms of the intrinsic cellular parameters, which are difficult to study biologically. We find, however, disagreement between novel biological data and our model in terms of the microenvironmental changes. We conclude that future theoretical models of stem-cell driven tumours must include specific feedback from the microenvironment onto the individual cellular behavior. Further, we identify several cell intrinsic parameters which govern loss of homeostasis into a state of uncontrolled growth.
As the importance of TICs becomes more evident as it pertains to aspects of tumour progression like heterogeneity [15], treatment resistance [16,17], recurrence [18] and metastasis [19], the need for generalizable therapeutic strategies based on conserved motifs in these cells grows. We therefore aim to understand how the phenotypic traits discussed earlier (asymmetric division rate, allowed rounds of transient amplification and lifespan of terminally differentiated cells) and microenvironmental changes (modelled as differences in oxygen supply) effect resultant tissue growth characteristics.
To this end, we present a minimal spatial, hybrid-discrete/ continuous mathematical model of a hierarchical SC-driven tissue architecture which we have used to explore the intrinsic, phenotypic, factors involved in the growth of TIC-driven tumours. We consider parameters that involve the rates of division of the cells involved in the hierarchical cascade as well as microenvironmental factors including space and competition between cell types for oxygen. We present results suggesting that there are discrete regimes in the intrinsic cellular parameter space which allow for disparate growth characteristics of the resulting tumours, specifically: TICs which form tumours that are unsustainable, TICs that are capable of forming only small colonies (spheres), and TICs that are capable of forming fully invasive tumours in silico, just as we see diversity in biological experiments ( Figure 2).

Results
A systematic parameter exploration of the three key parameters relating to vascularisation of the domain, symmetric vs. asymmetric division (a) and progenitor division potential (b) was performed. We also explored the parameter determining the lifespan of differentiated cells (c) and found that the only impact of longer lifespans is an increase in the amount of time before the simulations reach a steady state, but does not change the qualitative nature of the results. These results are summarised in Figure 3. Each of the three panels represents the results for a different degree of vascularisation (0.01, 0.05 and 0.1). A density of vascularisation of 0.05 would mean 12,500 oxygen sources in the domain. To determine the diffusion coefficient, we used the estimate of approximately 70 micrometers of effective oxygenation [20]. Each plot shows the total tissue size after 50,000 time steps as we change the proliferative potential of progenitor cells. Each of the lines shows a different rate of symmetric vs. asymmetric divisions. These results show that all these three parameters have a critical range where homeostasis is disrupted (tumourigenesis). Figure 4 shows examples of the typical results produced by this model. Although the proliferation rates of all the cells remain the same, due to space constraints and the differences in a, the population of TICs does not grow at the same rate as the non-stem population. Figure 4A shows an example of an unviable tissue (parameters: H~0:01, a~0:3, b~50 and c~1 day) where the vascularisation does not support the potential tissue size of that TIC, resulting in an area of hypoxia affecting the region that contains the TIC. That leads to the death of the stem cell and, eventually, the rest of the cells in the tissue. Figure 4B shows a case of slightly increased symmetric division, resulting in a dynamic homeostasis where cell birth and death is balanced so that tissue Figure 2. Differential phenotypes in cultures enriched for brain tumour initiating cells. Bright field images of CD133+ patient derived glioblastoma cell lines cultured in Neurobasal supplemented with EGF, FGF and B27, exhibiting striking phenotypic variability. These differences highlight the heterogeneity present even in a highly controlled static environment between cells that are putatively the same. doi:10.1371/journal.pcbi.1003433.g002 size remains relatively constant -which could represent the enigmatic dormant phase [9]. Finally, figure 4C shows an example where the system never achieves true homeostasis. In this case a is slightly higher when compared with the previous example, suggesting a critical value at which overgrowth occurs. Over time, the number of TICs increases, allowing for the 'tumour phenotype': unconstrained growth. Although this leads to areas of hypoxia, cells survive in the periphery of the blood vessels and keep growing until they take over the entire domain. A plot of cell number versus time for each of these three examples are plotted in figure 5.
Unsurprisingly, the higher the vascularisation of the domain the greater the tissue size it can support. Past a certain threshold, however, the difference becomes negligible and more remarkably, the qualitative dynamics are unchanged by any change in the microenvironment. The same effect is evident in the other two parameters, the rate of symmetric vs. asymmetric division (a) of TICs and the proliferative potential of TACs (b). Regardless of the vascularisation, disruption of homeostasis only occurs when the proliferative potential of TACs (b) is below a maximum value of about 15. For values of symmetric division (a) above 0.3, the values for b in which this overgrowth occurs becomes even more restrictive with a range of approximately 10-15.
Interestingly, we observed a conserved decrease in overall tissue size for the highest value of symmetric division, a~0:5, when the progenitor cells were allowed only 5 divisions (b~5). We believe   figure 4. The black trace, representing the unsustainable simulation, grows quickly though never expands its stem population and then outstrips the available oxygen and collapses. The blue trace, representing the homeostatic simulation, reaches a critical size and then maintains a steady birth-death balance. The red trace, representing the tumorigenic simulation, settles into an effectively linear trace on this log-log plot, suggesting power law growth. doi:10.1371/journal.pcbi.1003433.g005 this phenomenon represents a situation where the tissue is not able to grow to its potential as the stem cells themselves occupy too much space, and never allow the progenitors to contribute as much as they could to the overall population. This is a supposition however, and deserves closer study. These results are summarized in figure 3.
Of note as well: in no simulation did we observe the 'tumour phenotype' for a value of av0:3, suggesting something akin to a 'phenotypic tumour suppressor' function for this parameter. As observed biologically [10], this rate is highly susceptible to changes in microenvironment, suggesting an extension of this minimal model to include the microenvironmental factors measured in that study. How to incorporate the changes observed in that study into a mechanistic HCA model however, is not trivial, and we reserve it for a future extension of this work. Further, our current model exists only in two dimensions. While our quantiative parameters are based on experimentally derived values, the claims we make are largely qualitative abstractions, however, we stress that the specific quantitative descriptions of cell fates are likely not yet accurate and could change if this model was in three dimensions.

Discussion
In this paper we have presented a simple two dimensional computational model of the HM of a TIC-driven tissue. Our results show that there are distinct regions in parameter space (that directly correlate to the intrinsic TIC phenotype space) that encode vastly different behaviour in the tissue (or tumour) arising from the TIC in question. These parameters represent different TIC phenotypes, and therefore do not represent any specific genetic mutation. In this way, we hope to generalise the intrinsic alterations which a TIC could undergo much in the same way that the hallmarks of cancer have generalised non TIC-specific alterations [21]: our end goal being the identification of treatment strategies to target these phenotypes to slow or stop the progression of TIC-driven cancers.
Because of the difficulties in understanding TIC specific traits in vivo, the biological data to support these conclusions remains sparse. There have been some carefully undertaken in vitro experiments on single TICs in glioblastoma, a highly invasive and malignant brain tumour, which suggest that TIC specific division behaviour (a in our model) is variable and changes based on environmental cues [10]. Further work has shown that the other microenvironmental cues, such as acidity [14] and hypoxia [22][23][24][25][26][27][28] can also alter the prevalence of the stem phenotype by utilising functional markers of stemness, but the mechanism for this increase is, as of yet, imperfectly understood.
Our simulations do not show a significant TIC population dependence on vascular density (H), a surrogate for hypoxia, or a change in stem composition (see Supplementary Table S1), suggesting a flaw in the model. To rectify this, future iterations of this model should include direct feedback onto the cellular parameters from the microenvironment. We aim to parameterize this dependence by specific in vitro experiments designed to quantify this effect, rather than just elucidate its existence. Other future developments of this model should take into consideration the emerging body of work suggesting that the proportion of TICs within a tumour is directly affected by therapy and not just physiologic growth factor controls [29]. There is now evidence in several cancers to suggest that radiation increases the size of the TIC pool. Specifically, in breast cancer, it has been shown that radiation therapy induces non-stem cancer cells to de-differentiate into TICs [30]. Further, experimental studies have shown radiation increases the TIC pool in glioblastoma [31], which has often been attributed to radiation resistance associated with differences in cell survival [16]. A new study by Gao et al. [32], however, has shown in silico and in vitro that radiation can effect the symmetric to asymmetric division rate (our intrinsic parameter a), yielding further clues about the mechanism of this TIC pool expansion.
Dedifferentiation due to treatment related microenvironmental factors has not yet been considered in any spatial theoretical models. Dedifferentiation due to 'niche' specific factors was studied by Sottoriva et al. [33], who, using an agent based model, reported findings similar to ours: that the microenvironment made no significant change to the overall tumour growth dynamics. Beyond this single spatial study, the concept of SC dedifferentiation is gaining more and more attention in conceptual theoretical treatments [34] and has been modelled with a deterministic ordinary differential equation system for a well-mixed population of cells [35].
We, as well as others, find that the HM of tissue growth does not completely capture all the necessary dynamics that characterise cancer growth -but there is still a great deal of understanding to be gained from studying this formalism. To this end, we have performed a study of the factors related to TICs driving this dynamic and have identified several key factors which promote increased growth of the resultant tumour. Motivated by Hanahan and Weinberg [21], who have simplified the myriad (epi)genetic alterations which a tumour can undergo into the hallmarks of cancer, we seek to distill the traits of TICs in a similar way. Specifically, our model suggests that the number of allowed divisions of TACs exhibits bounds outside of which tumour growth is unsustainable. This finding has been corroborated independently by recent work from Morton and colleagues [36]. Further, there is a specific balance of symmetric to asymmetric division which keeps tumours from overgrowing; almost acting as a phenotypic tumour suppressor. Indeed, changes in this rate have been recently hypothesized to underlie the increasing stem pool in glioblastoma after irradiation [32], and could also hold a key to understanding tumour dormancy [9].
In summary, we have presented a minimal spatial Hybrid Cellular Automaton model of the HM of a TIC-driven tissue in which we have explored generalised TIC phenotypic traits and have identified several key cellular parameters which influence the overall tissue behaviour. While our model does capture a number of salient phenotypic characteristics of TICs that seem to be conserved, it fails to capture the recently observed changes in stem fraction secondary to microenvironmental perturbations. This is an indication that any computational model of a stem-hierarchical tissue, or tumour, built from this point on must not only include the physical microenvironment, but also feedback from the microenvironment onto the specific cellular parameters encoded in the HM.
Therefore, this endeavour has identified the crucial point that the microenvironment must effect the behaviour of the cells within the HM, and also several conserved phenotypic hallmarks, which could be the result of any number of (epi)genetic alterations or microenvironmental perturbations. By focussing on mechanisms important for the HM of stem-cell driven tumour growth, we are seeking to identify common phenotypes which could be targeted in a variety of solid tumours in which TICs promote tumour maintenance -thereby reducing the number of therapeutic targets to a more tractable set. Only with this sort of distillation of the biological complexity inherent to cancer initiation (and indeed progression) can we hope to make progress against this disease.

Materials and Methods
Our model is based on a hybrid, discrete-continuous cellular automaton model (HCA) of a hierarchically structured tissue. HCA models have been used to study cancer progression and evolutionary dynamics since they can integrate biological parameters and produce predictions affecting different spatial and time scales [15,33,[37][38][39][40]. As shown in figure 6C, cells are modelled in a discrete fashion on a 5006500 2-D lattice. This comprises approximately 1cm 2 where we assume a cell diameter of 20 micrometers [41]. The domain has periodic boundary conditions but the simulations are stopped when a cell reaches one of the boundaries. Every time step, cells are iterated in a random fashion as to avoid any bias in the way that cells are chosen. Figure 6A shows that, although all cells are assumed to have the same size and shape, they can only be one of three different phenotypes: TICs capable of infinite divisions, TACs which are capable of division into two daughters for a certain number (b) of generations, and TDs which cannot divide but live and consume nutrients for a specified lifetime (c). Modes of division for TICs include asymmetric division (with probability 1{a), which is division into one TIC daughter and one TAC daughter and symmetric division, which is division into two TIC daughters (probability a).
The continuous portion of this model is made of up the distribution and consumption of nutrients (in this case modelled only as oxygen). Vessels, which are modelled as point sources and take up one lattice point (V i,j in Equation 1), are placed randomly throughout the grid at the intiation of a given simulation, in a specified density (H). Each of these vessels supplies oxygen at a constant rate (l) which then diffuses into the surrounding tissue. The diffusion speed/distance is described by Equation 1, where O(x,y,t) is the concentration of oxygen at a given time (t), and place (x,y), D O is the diffusion coefficient of oxygen, l is the rate of oxygen production from a blood vessel, m s , m p , and m t are the rates at which TIC, TAC and TD cells consume oxygen. The difference in time scales that govern the diffusion of nutrients and that at which cells operate is managed by updating the continuous part of the model 100 times per time step. During each update the oxygen tension in a given grid point is updated with the values of the surrounding cells using a von Neumann neighbourhood modulated by the diffusion rate (D O ).
LO(x,y,t) Lt~D O + 2 O(x,y,t)zlV x,y {m S S x,y {m P P x,y {m T T x,y ð1Þ Any simulation performed by this model can be characterised by the parameters found in figure 7. The most relevant parameters for the question we are trying to address are the following:  In each case, as can be seen in figure 6, a simulation is seeded with one TIC with a given set of intrinsic parameters (a, b, c) governing its and its progenys behaviour, which is placed in the centre of the computational domain. The domain is initialised with as many randomly placed oxygen source points (vasculature) as described by the vascular density parameter (H).

Supporting Information
Table S1 Raw stem and total cell numbers from several runs of the CA with varying parameter combinations. (XLSX)