Cancer Stem Cell Tumor Model Reveals Invasive Morphology and Increased Phenotypical Heterogeneity

Malignancies arise after sequential accumulations of mutations in oncogenes and tumor suppressor genes (1). During the process of malignant transition, fundamental regulatory mechanisms are lost. The cancerous cell population has unlimited growth potential, evades the immune system and apoptosis, and often acquires the ability to breach tissue boundaries and expand into foreign environments (1). The cancer stem cell (CSC) concept sheds new light on all of these features (2, 3). In this study, we investigate the influence of a hierarchical organization of malignant clones on tumor growth, evolution, invasion, and morphology using a multi-scale cellular automaton–based computer model (4).

Malignancies are highly heterogeneous tissues containing largely diverse cancer cell populations as well as other cells such as fibroblasts and macrophages (1). The dominant genetic view of malignancies explains the heterogeneity in cancer cells with the presence of genetically diverse clones emerging from the continuous acquiring of additional genetic lesions by cancer cells. Such clones compete for resources, resulting in a highly dynamic process known as Tumor Darwinism (5). In this article, we refer to this view of malignancies as the classical model. Although this model greatly contributes to our understanding of malignancies, recent experimental evidence suggests an additional layer of complexity. The heterogeneity present in tumors could, in part, be the result of the diversity in differentiation grade of genetically identical cells (3, 6). For example, in glioblastoma multiforme, cells with an immature phenotype expressing the cell surface marker AC133 are the cells that fuel tumor growth and have the exclusive capacities to self-renew, differentiate, and transplant the malignancy into severe combined immunodeficiency mice (7). Self-renewal and spin-off of differentiated cells are features shared with normal stem cells, and therefore, cells with such features in malignancies are defined as CSCs (6, 8).

Tumor growth is generally accepted to be the result of several highly complex interacting processes. Fundamental cellular characteristics such as genetic and epigenetic features influence signal transduction route activities that in turn control cellular functions such as mitosis, apoptosis, and cell migration. In addition, environmental factors including nutrients and growth factor concentrations interplay with these processes. To study the emergent properties of such systems regarding proliferation speed, infiltrative growth, and phenotypical evolution of cancer, advanced mathematical models have been developed (9–11). Here, we apply computational modeling techniques to investigate the consequences of hierarchically organized cancer cell populations on solid tumor growth dynamics and progression. We describe that implementing the developing concept of CSCs in a mathematical tumor growth model directly results in an invasive morphology. Moreover, we found that hierarchical organized malignant clones have highly altered evolutionary dynamics. Most strikingly, the CSC organization promotes phenotypical heterogeneity, a feature that could have immediate consequences for therapy resistance.

We developed a hybrid tumor growth model based on cellular automata (4) and partial differential equations. We refer to this model as the SCA model. In the SCA model, the individual cancer cell is the fundamental unit of the tumor, we simulate its proliferation, metabolism, migration, stemness, and differentiation (see Quick Guide and Supplementary Materials and Methods).

Implementing the CSC model of malignancies means to simulate cancer cells with different replicative potential within the tumor. For simplicity, we assume that in our model there are only two types of cells: CSCs and DCCs. CSCs possess unlimited replicative potential and could either generate new CSCs (with a probability P_S) or DCCs. DCCs can divide for a maximum of H generations before stopping to proliferate irreversibly. This method yields a hierarchy with CSCs at the top and DCCs at the bottom. We simulate the classical model of malignancies, in which all cells possess tumor growth–promoting capacities, by simply setting P_S = 1. In such a situation, all cells possess stem cell characteristics. Using this method, we have an intuitive way to compare the flat, classical tumor model with the hierarchically organized CSC model.

Computational modeling allows the exploration of highly complex systems, such as tumor growth. In this study, we have used a computational tumor growth model to test the consequences of hierarchical organized clones on a vast range of areas of tumor biology. These include tumor growth dynamics and morphology, but also tumor evolution and therapeutic effects. The applied computational modeling technique allows us to get more insight into the underlying dynamics of these aspects of cancer growth and progression that would be impossible in a conventional experimental biological setting.

First, we investigated how tumor growth dynamics are altered upon varying the CSC fraction. In the SCA model, this corresponds to changing the variable P_S. For high values of P_S, we expect to model the classical interpretation of tumors because all cell divisions are symmetrical and all cells are therefore clonogenic. In contrast, low values of P_S represent the CSC model. We simulated the growth of tumors with P_S = 1, P_S = 0.1, and P_S = 0.03 with the variables described in Table 1.

Figure 1A shows the fraction of CSCs on the total amount of tumor cells for different values of P_S. The selected P_S values correspond to CSCs populations comprising roughly 100%, 1%, or 0.1% of the total tumor volume and therefore cover mainly the CSC fractions observed in a variety of solid malignancies (2). In the CSC model (P_S = 0.1 and P_S = 0.03), small (early) lesions have relatively high fractions of CSCs, although this number decreases and tends to stabilize when the tumor progresses. This observation is supported by in vivo studies that find increased numbers of CSC marker–bearing cells in micrometastases compared with larger tumors (21). This indicates that even with fixed self-renewal rates (P_S), this phenomenon is intrinsic to lesions initiated by a single CSC although environmental factors influencing self-renewal frequencies are also likely to contribute. Tumor growth curves for various P_S values all display the classical Gompertzian-like growth kinetics. However, as expected with equal cell cycle durations, the self-renewal rates of the stem cell fraction influences proliferation rate greatly, hence, low self-renewal rates (small P_S) correspond to slow tumor growth (Fig. 1B). Interestingly, the decrease in accumulation of tumor volume is accompanied by a stabilization of the fraction of CSCs, suggesting an intimate relationship between these two processes.

Spatially, all experiments display a three-layered structure consisting of an external proliferative area, an inner senescent layer, and a necrotic core. However, tumor morphologies for different P_S values are remarkably dissimilar (Fig. 1D). For P_S = 1, in which there is no hierarchy, a symmetrical, sphere-like tumor morphology is generated that closely resembles early tumor growth models (22, 23). In contrast, the shape of the tumors generated with low P_S values is highly irregular (especially P_S = 0.03). CSC-driven tumors yield highly invasive morphologies with fingering fronts and clusters of cancer cells beyond the tumor margin, driven by the mobility and the exclusive proliferative properties of CSCs (Movies S1 and S2).

From Fig. 1C, it is evident how hierarchically organized tumors (P_S = 0.03 and P_S = 0.1) generate a higher degree of invasiveness (see Supplementary Materials and Methods), compared with tumors in the classical model (P_S = 1). It is important to note that the intrinsic properties of the cells in the classical tumor and the stem cells in the CSC-driven tumor are completely identical.

Magnification of a tumor border in a CSC-fueled tumor growth model (P_S = 0.03) shows how CSCs migrate beyond the margins of the tumor mass. CSCs colonize the surrounding tissue and expand locally, forming small satellites that grow back into, and are engulfed by, the main tumor mass (Fig. 2A). This is accompanied by the observation that the invasive behavior of individual experiments that make up Fig. 1C occurs in a wave-like fashion (Supplementary Fig. S1). These results are paralleled by recent findings in a different model system in which high migration levels in a small subset of cells give rise to small proliferating extratumoral lesions and therefore tumors are “conglomerates of self-metastasis” as the authors propose (24).

In an endeavor to validate these observations from our computational model, we investigated human tumor specimens. Close examination of the histology of different highly diverse human malignancies displayed a relatively confined large tumor mass with clearly detached tumor cells forming small lesions in the surrounding normal tissue (Fig. 2B). This exemplifies that human tumor histology contains indications of a stepwise infiltration and colonization of the surrounding tissue as our model predicts would follow from a hierarchical organized malignancy. Next, we attempted to determine the relationship between CSC fraction and invasive properties, as predicted by the model, using in vitro cell culture (Fig. 2C; Supplementary Figs. S2–S5 and Table S1). We determined the clonogenicity of a set of various cell lines (n = 24) as a surrogate for their CSC fraction. Additionally, we quantified the invasive properties of clonally derived structures with our measure of irregular morphology both on adherent plates and in Matrigel for all these lines. A significant inverse relationship between clonogenic fraction and the invasive properties of these lines exists. This implicates that tumor structures driven by a small fraction of clonogenic cells tend to generate a more irregular and invasive morphology, a finding that corroborates the predictions of our model. Combined, we take this as evidence that the SCA model based on the CSC concept closely resembles tumor growth patterns and morphology.

To show that our results are not limited to two dimensions only, we expand the SCA model to three dimensions (see Supplementary Materials and Methods). The three-dimensional SCA model reproduces the invasive tumor morphology induced by the hierarchical organization, yielding three-dimensional fingering tumor fronts and clusters of cancer cells beyond the borders of the main tumor mass (Fig. 3A; Movie S3). On the contrary, again a classical model with equal volume does not display any apparent invasive pattern and instead exhibits a spherical morphology (Fig. 3B; Movie S4). We therefore propose that the CSC model, in contrast to the flat model, more faithfully reproduces the human tumor morphology for both the two-dimensional and the three-dimensional implementations.

The dominant model evoked to explain the advancement of malignancies is clonal evolution, as was first proposed by Nowell (5). The term effective population size is used in population genetics to indicate the fraction of total individuals in a population that effectively contributes to the next generation and are therefore evolutionarily relevant. Hence, the effective population size in a CSC-driven malignancy is smaller than in the classical model because only mutations in the CSC compartment contribute to the evolutionary process (2). To implement clonal evolution in our model, we assume that, at each symmetrical division a CSC has a probability P_Mut to acquire a genetic hit and generate a daughter cell with a different phenotype selected from a randomly generated pool of 30 phenotypes (Table S2 and S3; see Supplementary Materials and Methods). Under an equal mutation rate (P_Mut = 0.1), the CSC model exhibits a slower acquisition of new phenotypes due to its smaller effective population size (Fig. 4A), but also shows a radically different selection process. Strikingly, although the rate of emergence of new phenotypes is much lower than observed in the classical model, a wide range of newly acquired phenotypes expands and contributes to the malignancy (Supplementary Fig. S6A and B).

To compare the phenotypical selection more closely, we synchronize the pace at which new traits emerge. Evolutionarily, a CSC model with P_S = 0.03 and P_Mut = 0.5 corresponds to a classical model with P_S = 1 and P_Mut = 0.001 (Fig. 4B). Intriguingly, despite this adaptation of the mutation pace, the two models differ substantially in the way they exercise clonal selection. From a representative example in Fig. 4C, it is evident that the CSC model allows for a much higher phenotypical heterogeneity whereas the classical model seems to select for a small number of aggressive clones (average result in Supplementary Fig. S6C). Quantification of the heterogeneity clearly underscores this (Supplementary Fig. S7A and Supplementary Materials and Methods). It is important to realize that the selective pressure from the environment is equal in both models. This, combined with the equal rate of new phenotype occurrence, suggests that the geometric properties of the CSC model promote heterogeneity. We argue that the intrinsic properties of the CSC model might propel an alternative process to natural selection, referred to as genetic drift (25). In populations with small effective population sizes, sampling errors are frequent and might allow for the expansion of clones with no clear survival benefit. Interestingly, the invasive properties of the CSC model might fuel such a sampling error promoting mechanism. The phenomenon we observe at the tumor margins during invasion: a CSC migrating out of the tumor initiating a small satellite lesion, resembles the Founder Effect (26). The Founder Effect occurs when a new population is established by a low number of individuals from a larger population, this process is often accompanied by a loss of phenotypical variation in the new population due to sampling errors. In such a scenario, the new generation can differ substantially from the previous and is not in direct competition for nutrients and space. This implies that the pattern of tumor growth in the CSC model stimulates genetic drift at the infiltrative edges of the tumor and therefore promotes phenotypical heterogeneity.

After having established the effects of CSC-driven tumor growth on the invasion and evolution of malignancies, we now explore the crucial topic of therapeutic intervention and tumor relapse. CSCs have been suggested to be more resistant to therapeutic interventions such as chemotherapy or irradiation compared with their differentiated counterparts (27, 28). Also tumors that relapse after seemingly successful therapy are believed to regrow from the CSC that survived the therapeutic regimen (27). Here, we investigate the dynamics associated with therapeutic interventions that are either selective for non-CSCs or equally efficient against both cell types. We find that the morphology and growth kinetics of relapses for both types of therapeutic interventions are very much different (Fig. 5; see Supplementary Materials and Methods). Regrowth after therapy that specifically targets DCCs is accompanied by enhanced invasive growth patterns whereas relapsing tumors after stochastic tumor cell killing are similar to the malignancy before treatment (Fig. 5A). Simultaneously, in case CSCs are resistant to therapy, the pace at which the malignancy relapses is greatly enhanced due to the relatively high fraction of CSCs directly following therapy (Fig. 5B). Also, the invasiveness of the recurrent tumors is markedly increased following intervention which is not effective against CSCs (Fig. 5C). These findings are in line with a range of clinical observations describing increased growth speed and enhanced invasion in the relapsing malignancy that are mostly attributed to the selection of more aggressive clones by the drug (29–31). However, our data now indicate that these observations could be partially explained by the failure of conventional therapies to eradicate the CSC compartment and the subsequent relapse dynamics in CSC-driven tumors. Evaluation of evolutionary dynamics during relapse after both types of intervention revealed significant differences as well. Following therapy which is ineffective against CSCs, relapsing tumors display a marked increase in heterogeneity, whereas therapy that does target CSCs results in a dramatic decrease of heterogeneity (Supplementary Fig. S8). This latter scenario is related to the fact that relapses are very much different compared with the primary malignancy with respect to the clonal lineages that contribute to the relapse of the tumor (Fig. 5D). Combined, this implicates that applying therapy that is ineffectively targeting the CSC population is not only unsuccessful in curing the patient but also promotes malignant features including rapid expansion, increased invasion, and stimulates heterogeneity directly after therapy.

We have applied a computational tumor growth model to investigate the effects of hierarchical organization on a range of areas in tumor biology. The obtained results provide novel insights into the influence of CSC-driven tumor growth on tumor morphology and invasion. Importantly, we have been able to partially validate these findings in a biological experimental setting. We have also established a significant effect of hierarchical organized malignant clones on tumor evolution, and consequently, on tumor progression and relapse after therapy. Our results stress the need to develop therapeutic interventions that efficiently target the CSC compartment because drugs that fail to do so are not only unable to cure the patient—but even enhance the malignant properties of tumors. Therefore, we conclude that computational oncology has the potential to greatly enhance our understanding of fundamental cancer biology and improve future therapeutics.

The capacity to display infiltrative growth is a pathognomic feature of malignant cells (1). Crossing tissue boundaries is the first step in the process of metastasis and is therefore of special interest. Surprisingly, this aspect of tumor growth has rarely been investigated in mathematical models. However, recent important work of Anderson and colleagues (11) and Bearer and colleagues (32) illuminate how selective pressure, applied by the microenvironment, in combination with tumor cell heterogeneity, could result in an invasive morphology. Here, we report how, in the SCA model, invasion emerges from the hierarchical organization of malignant clones and need not be driven by external interaction with a heterogeneous microenvironment. Hence, invasive morphology is an intrinsic feature of CSC-driven malignancies and has to be considered when studying tumor invasion.

Simultaneously, high levels of clonal heterogeneity within malignancies are associated with rapid progression of the disease (33), poorer survival (34), and occurrence of therapy resistance (35). Understanding this phenomenon in more detail might improve the design of more effective treatments and chemopreventives. Our results show how CSC-driven tumor growth stimulates heterogeneity in malignancies. This is related to the reduced effective population size in these malignancies and the typical growth pattern that results in the segregation of clones by formation of small invasive lesions. These satellite structures are not in direct competition with the dominant clone, allowing suboptimal clones to expand. We speculate that interfering with these dynamics, for example by inhibiting migration, and thereby segregation, might be a potent means by which reduced heterogeneity accompanied by the slower progression of the disease and better response to therapy can be achieved either in a chemopreventive or therapeutic setting.

Additionally, our results underlie the need to develop therapeutic modalities that successfully eradicate the CSC compartment. We show how therapeutic interventions that do not target CSCs efficiently are not only ineffective as a cure, but might even contribute to an increase in malignant properties of the relapsing tumor. We find, for example, enhanced heterogeneity and invasion in the relapsing tumors that are treated with interventions that are only directed against the differentiated cells. This is not due to the selection of the most malignant clones present in the cancer tissue but are intrinsic to the dynamics of regrowth in CSC-driven malignancies.

The model presented provides a range of predictions and implications regarding CSC-fuelled tumor growth that can be tested and exploited to investigate the properties of hierarchical organized malignancies. Here, we report that the clonogenicity of cell lines is connected with the invasive properties of these lines in vitro as suggested by our model. In future research, it would be interesting to expand this finding to human tumor specimens or in established CSC lines in which the CSC fraction can be manipulated. Additionally, our model illuminates how tumor features, such as clonal heterogeneity, are related to the CSC fraction (Supplementary Fig. S7B). Therefore, we speculated that careful analysis of these tumor properties might lead to optimized techniques that establish the CSC fraction in malignancies without the need for transplantation studies which are currently highly disputed (36, 37).

The current version of the SCA model clearly shows the dynamics of CSC-driven tumor growth and its consequences for tumor morphology and evolution. However, further research efforts will undoubtedly lead to increased insight into the nature of the hierarchical organization of tumor cells. If so, the SCA model can be easily adapted to implement potential new information and subsequently come to even more accurate description of tumor growth. For instance, we have assumed several variables such as CSC self-renewal frequency, mutation rate, and migration speed to be constant. However, in in vivo tumors, these variables might be more dynamic both throughout tumor progression and as a consequence of microenvironmental interactions. For example, CSCs have been suggested to reside in a so-called “CSC niche” that is formed by secreted factors and direct cell-cell interactions (38, 39). This niche is believed to protect CSC from differentiation and to promote CSC self-renewal and therefore introduces dynamic P_S values. Moreover, physical circumstances such as oxygen concentration might have differential effects on CSCs and DCCs, and therefore, might play a role in dynamic CSC functions (40). In addition, mutation rates are potentially subject to changes during tumor development (41) and intriguingly, might even be influenced by environmental conditions, such as that reported for bacteria (42). Currently, it remains elusive as to what extent this latter phenomenon, which is referred to as “adaptive mutation,” is involved in human malignancies (43). All these considerations justify further investigations, but will be dependent on thorough experimental examination of the variables involved.

Importantly, from the current formulation of the SCA model, it is indisputable that hierarchical organization of malignancies significantly contributes to the invasive morphology and increased heterogeneity of tumors and is therefore a crucial issue for better understanding tumor biology and to improve current anticancer treatments.

No potential conflicts of interest were disclosed.

We thank M.R. Sprick, D.J. Richel, and F. de Sousa Mello for useful comments and R. Belleman and J.H. de Jong for technical assistance.

Grant Support: Academic Medical Center (J.J.C. Verhoeff, J.P. Medema, and L. Vermeulen), ZonMW VICI program (J.P. Medema), and the Faculty of Science of the University of Amsterdam (A. Sottoriva, L. Naumov, and P.M.A. Sloot).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.