Selection-driven tumor evolution with public goods leads to patterns of clonal expansion consistent with neutral growth.

Cancers are the result of eco-evolutionary processes fueled by heritable phenotypic diversification and driven by environmentally dependent selection. Space represents a key growth-limiting ecological resource, the ability to explore this resource is likely under strong selection. Using agent-based modeling, we explored the interplay between phenotypic strategies centered on gaining access to new space through cell-extrinsic degradation of extracellular matrix barriers and the exploitation of this resource through maximizing cell proliferation. While cell proliferation is a cell-intrinsic property, newly accessed space represents a public good, which can benefit both producers and non-producers. We found that this interplay results in ecological succession, enabling emergence of large, heterogeneous, and highly proliferative populations. Even though in our simulations both remodeling and proliferation strategies were under strong positive selection, their interplay led to sub-clonal architecture that could be interpreted as evidence for neutral evolution, warranting cautious interpretation of inferences from sequencing of cancer genomes.

Introduction

Cancer starts when a somatic cell breaks away from the rules of homeostatic cooperation imposed by multicellularity and starts acting as an independent evolutionary entity. As this cell divides, its progeny (clone) acquires genetic and epigenetic alterations generating novel phenotypic variants (subclones). This subclonal diversification in the context of competition for growth-limiting resources between cancer cells enables Darwinian evolution, which can lead to the emergence of key tumor phenotypes, defined by hallmarks of cancers, culminating in lethal malignancies (Axelrod et al., 2006; Greaves and Maley, 2012; Hanahan and Weinberg, 2000, 2011; Korolev et al., 2014; Nowell, 1976). Consistently, genetic and phenotypic intratumor heterogeneity (ITH) correlates with poor prognosis and contributes to treatment failure (Gerlinger et al., 2012; Swanton, 2012). Understanding the basic rules that govern the somatic evolution and how ITH is maintained in the face of subclonal competition is thus key to a more effective tumor prevention and therapies (Basanta and Anderson, 2013; Sottoriva et al., 2015).

Mathematical modeling has been used to understand the behavior of complex dynamic systems, such as changes in cell populations in space and time during the emergence and progression of cancers, which are difficult to intuit without the help from rigorous theoretical frameworks. Properly parameterized and experimentally tested mathematical models can provide a deeper understanding of these complex processes and also serve as predictive tools for therapy optimization. Ordinary differential equations, the most commonly used tools to model somatic evolution, typically assume well-mixed populations (Altrock et al., 2015; Gatenby and Gillies, 2008). This assumption fails to account for heterogeneity in space and tissue architecture as a key growth limiting resource, which has a profound impact on evolutionary dynamics. Agent-based models (ABMs), which can consider cells as individual agents rather than populations, can be used to account for space and spatial structures. The importance of space and the applicability of ABMs has been explored showing that an increase in available space results in higher ITH and can affect evolutionary mode (Chkhaidze et al., 2019; Noble et al., 2019; West et al., 2019).

Epithelial tissues, which give rise to the majority of human cancers, are organized as layered structures, constrained by basal membranes and the extracellular matrix (ECM) (Bissell et al., 2002; Bissell and Radisky, 2001). Carcinogenesis involves breaking from the growth limitations imposed upon by epithelial layers and gaining access to the new space through the degradation and remodeling of the ECM (Cawston and Wilson, 2006). From an ecological perspective, accessing this new space constitutes niche engineering; thus, we will refer to the phenotypes of cancer cells which possess the ability to remodel the ECM as “engineers” (Lloyd et al., 2016; Myers et al., 2020). ECM degradation is achieved through the secretion of matrix-degrading enzymes, which require a sufficiently high local concentration to achieve their effects. Since the newly accessed space can benefit not only engineers but also non-engineer cells, the space, unlocked by engineers, represents a “common good” resource (Kagel and Roth, 1995). Optimization for use of space can be achieved by the increased consumption of this resource through elevated proliferation rates; thus, we will refer to the phenotype that optimizes space consumption as “proliferators”. To understand eco-evolutionary dynamics resulting from selection for optimal phenotypic strategies toward creating (engineering) and consuming (proliferating) space, we consider an abstract on lattice grid, where the domain is divided up into multiple subdomains by degradable ECM barriers.

Using a cellular automaton to implement our ABM, we examined the eco-evolutionary dynamics emerging under selection pressures that act to optimize mutually exclusive engineering and proliferating phenotypic strategies within domains with different spatial organization (Gatenbee et al., 2019a; Margolus and Toffoli, 1987; Poleszczuk and Enderling, 2014). We found that gaining access to new space through non-cell autonomous effects of engineers can dramatically increase ITH and lead to patterns which, in the context of ecology, can be described as ecological succession. Surprisingly, we found that this selection-driven dynamic could lead, under sampling resolutions commonly used in cancer genomics, to subclonal diversification patterns that could be misinterpreted as evidence for neutral evolution.

Results

Model description

We initiate simulations with 25 cells with the basic phenotype seeded in a subdomain, surrounded by ECM barriers. The subdomain itself is located in the center of the large domain, composed of the grid of subdomains. Cell growth is constrained to the subdomain, unless the ECM separating the population from the neighboring subdomain is degraded. ECM barriers can be degraded through action of engineer cells if a sufficiently high concentration of ECM-degrading enzymes is produced – which is a function of both the number of engineers, as well as their efficiency (Figure 1A). When a cell divides, it can acquire a “driver” mutation that permanently increases either expression of ECM-degrading enzymes (for the engineer phenotype) or proliferation rate (for proliferator phenotype). After each division, a check is made to remove any ECM which has met the removal criteria. A cell can die through apoptosis (with a rate of 15% each timestep), freeing up space that can be explored by its neighbors. There is phenotype plasticity in the model, so engineers can become proliferators (at a rate of 0.2%) and vice versa. Although these engineering and proliferative phenotypes are mutually exclusive, the mutational history of each cell enables them to start from the previous rates of engineering and proliferation if they switch back (Figure 1B). Simulations are initialized within an 601x601 domain, separated into smaller sub-domains, ranging from 25x25 to 301x301 (shown in Figure 1C), and run for 2500 timesteps, which allows for simulations to run for sufficiently long time so that the relevant evolutionary dynamics can be examined.

Figure 1

Model description

(A) A depiction engineers degrading an ECM barrier in the model. For the ECM to be degraded, a sufficient concentration (~45% enzyme per lattice unit) of matrix-degrading enzyme is required near the ECM. The top images show the concentration of enzyme and the bottom images show the population of engineers producing them (with phenotypic color map below).

(B) A depiction of how driver mutations interact with phenotypic switching in the model. The top plot shows the proliferation driver mutation space, and the bottom plot shows the engineering driver mutation space. For all basic cells, their first driver mutation is randomly selected, which is the only time a mutation determines the cell's phenotype (orange box). Afterward, driver mutations will increase the cell's phenotype at the time of the mutation. However, phenotypic switches occur at the end of a time frame. The change boxed in blue shows the cell acquiring a proliferative driver (as it was the cell's phenotype at the time of division), yet it stochastically switched phenotypes at the end of the timestep. Similarly, the cell acquires an engineering mutation in t = 3 and again stochastically switches the phenotype at the end of the timestep (boxed in pink). The cell retains a memory of its mutational history, so the fitness does not reset upon switching phenotypes, yet only one phenotypic characteristic can be expressed at a time.

(C) Images taken from simulations of the four subdomain sizes in the model, showing the two types of characterized simulations. All parameters are the same except subdomain size. The left column shows simulations, where engineering failed, and the right column shows simulations where engineering succeeded.

Impact of environmental engineering on phenotypic and mutational diversity

Regardless of the subdomain size, we observed two scenarios characterized by engineering success. In some simulations, the engineering phenotype was driven to extinction, which left the population confined within the initial subdomain (failed engineering). In others, the population was capable of accessing the majority or all of the subdomains of the larger domain (Figure 1C). As a control, we examined the impact of domain size on phenotypic and clonal heterogeneity in the complete absence of engineering phenotype (Figure S1A). In the absence of successful engineering, phenotypic diversity peaked early and then decreased as more effective proliferators were selected. Increase in the subdomain size lead to an increase in the peak phenotypic diversity. In the presence of successful engineering, phenotypic heterogeneity peaked at later time points but at significantly higher levels. Similar to the outcome in the absence of engineering, the phenotypic diversity decreased after reaching its peak, indicating convergent evolution toward highest proliferation rates (Figures 2A and S1A). Clonal heterogeneity increased directly with subdomain size, reflecting an increase in population size (Figures 2B and S1B). After peaking, upon population reaching carrying capacity within spatial constrains, clonal heterogeneity remained high, with the effect more easily observed in larger subdomains. Successful engineering led to a delay in reaching maximal clonal heterogeneity but enabled populations to reach much higher levels of clonal heterogeneity. Notably, engineering limited reduction in clonal ITH, observed shortly after the peak. Therefore, environmental engineering resulting in production of public goods (space) enabled tumor populations, not only to gain access to space and achieve larger population size but also to gain increased levels of phenotypic and clonal heterogeneity.

Figure 2

Impact of successful engineering on intra-tumor htheterogeneity

Temporal changes in genetic (A) and phenotypic (B) heterogeneity are captured by the Shannon Index. Green and black lines represent successful and failed engineering, respectively. Each plot depicts the outcomes of 80 simulations. NS, ∗, ∗∗, ∗∗∗, and ∗∗∗∗ represent p values of Kolmogorov-Smirnov tests at the indicated time points: NS indicates p > 0.05, ∗ indicates p < 0.05, ∗∗ indicates p < 0.01, ∗∗∗ indicates p < 0.001, and ∗∗∗∗ indicates p < 0.0001.

Successful engineering was associated with noticeable patterns of phenotypic succession. After breaking the initial ECM barrier, engineers exploited the newly accessed space, growing outward the next barriers. As engineers continued to degrade the ECM and grow outward, proliferative cells began to outcompete vestigial engineers near the core of the tumor (Figure 3A, Video S1). Once all of the ECM barriers in the domain were removed, engineers were succeeded by proliferators, both from the expansion of pre-existing proliferative sub-clones and mutational conversion of engineers to proliferators (Figures 3B and 3C). Notably, this phenotypic succession was not accompanied by a clonal one, as while phenotypic strategy converged toward proliferative phenotype, clonal heterogeneity was maintained at near peak levels. This result implies high levels of convergent evolution, where identical phenotypic solutions are reached by different subclonal lineages.

Figure 3

Patterns of ecological succession in the model

(A) Images from the visualization of a simulation outcome showing engineers gaining access to new space and being succeeded by proliferative cells over time.

(B) Frequencies of basic, proliferative, and engineering cells over time.

(C) Populations of proliferative cells, engineering cells, and total tumor over time. Engineers flourish when they break out initially, and over time, proliferative cells can grow into the space created by engineers.

Video S1. The video shows what the model visually looks like when executed

For simplicity, the model’s current state is output as a .png image every 10 timesteps. These images have been compiled together into a short video to show tumor growth over time, and how engineers degrade the ECM.

Selection-driven dynamics can lead to clonal architectures consistent with neutral evolution

Since we recorded mutational history for each of the cell within simulation, these data enabled reconstruction of “true” (within simulations) clonal architecture, in contrast to inferred clonal architecture obtained from analyses of bulk genome sequencing. We visualized clonal architecture using Muller plots (Gatenbee et al., 2019b). In this visualization, subclones with sub-threshold frequency (the majority of subclones in our simulations) are “invisible”, instead being grouped with the parental (sub)clone. We chose to visualize changes in clonal architecture using 10% clonal resolution threshold (Figure 4A), corresponding to a common 20x genome sequencing depth (assuming near diploid genome), and 1% clonal resolution threshold, reflecting higher resolution analyses (Figure 4B), and 0.1% threshold (Figure S2) – which provides a more accurate representation of the “true” clonal architecture.

Figure 4

Muller plot visualization of the spatiotemporal clonal dynamics

Two representative simulations for each subdomain size are shown, one where engineering failed and one where engineering succeeded. The x-axis on each plot is time, and the y-axis is population size. Colors indicate time of clone's emergence, with earlier clones presented by dark blue and later clones presented by warmer colors. Only clones with above-threshold frequency are plotted.

(A) Visualization with a frequency threshold of 10%.

(B) Visualization with a frequency threshold of 1%. Visualization with a frequency threshold of 0.1% is shown in Figure S2. Plots at different resolution levels represent the same individual simulations.

In the absence of successful engineering, in the smallest subdomain size, we observed patterns of clear clonal succession (Figure 4). Increase in the subdomain size resulted in tumors with more complex clonal architecture but still with clear evidence of clonal expansions. Successful engineering dramatically increased clonal heterogeneity, leading to lower rates of expansion of individual sub-clones, indicative of increased clonal interference. Noticeably, increase of the analysis resolution revealed more complex clonal patterns (Figures 4 and S2).

Presence/absence of strong clonal expansions within the Muller plot visualization can be intuitively interpreted as evidence of selection-driven/neutral evolution. However, as such an interpretation might be misleading, we asked whether the apparent differences in patterns of clonal expansions observed at different resolution levels would be interpreted differently when analyzed using quantitative metrics of neutrality. To this end, we examined the relation between the cumulative number of mutations as a function of the inverse frequency of the mutations. This approach has been used to discriminate between neutral (no fitness differences between subclonal lineages) and selection-driven evolution in inferences of evolution modes from bulk sequencing data, where a linear fit between the inverse of the variant allele frequencies (VAFs) of the detected mutations with the cumulative mutation function is interpreted as evidence of evolution neutrality (Williams et al, 2016, 2018). Since we had VAF data recorded from our simulations, we analyzed the goodness of linearity fit of our simulation data at different sub-domain sizes and resolution levels of 1%, 0.1%, and 0.01%. As a reference point, we generated VAF data from simulations that followed neutral dynamics (no changes in proliferation probability), under a scenario of no ECM barriers (lack of subdomains). Strikingly, we found a better linear fit of data from our selection-driven simulations, compared to the neutral control, at both 1% and 0.1% resolution, with the trend reversed at 0.01% (Figure 5A). The impact of resolution level on goodness of fit was not limited to our simulations, as we have also observed marked differences in the goodness of linearity fit within the artificial data set provided with neutralitytestr (Williams et al, 2016, 2018) (Figure S3). Interestingly, comparison of the goodness of linearity fit between simulations with successful and failed engineering revealed significantly higher fit in simulations with successful engineering across all of the resolution levels, with the differences being more pronounced at lower subdomain size (Figure 5B). These results warrant caution in inferring mode of evolution from genomic sequencing data.

Figure 5

Consistency of clonal architecture with metrics of evolution neutrality

Goodness of fit to a linear regression between the cumulative mutation function and the inverse VAFs (M(f) ∝ 1/f) is used as a metrics of evolution neutrality.

(A) Goodness of feet for the in-silico growth simulations from Figure 2 with the indicated domain size and simulations with neutral growth dynamics. Statistical comparisons are made between simulations with successful engineering and the neutral simulations.

(B) Impact of engineering success on the goodness of fit. Statistical comparisons are made between the success and failure (indicated by + and – respectively) of each subdomain size. NS, ∗, ∗∗, ∗∗∗, and ∗∗∗∗ represent p values of Kolmogorov-Smirnov between two distributions: NS indicates p > 0.05, ∗ indicates p < 0.05, ∗∗ indicates p < 0.01, ∗∗∗ indicates p < 0.001, and ∗∗∗∗ indicates p < 0.0001.

Discussion

Our study examined the eco-evolutionary dynamics of tumor cell populations using an ABM where cells evolve under explicit pre-defined rules, specifying the impact of mutations and epigenetic switches on phenotypes that impact Darwinian competition for a limited ecological resource – space. In addition to the commonly used consideration of evolutionary fitness as a cell-intrinsic proliferation probability, our simulations included a second evolutionary strategy, which specified the ability of cells to unlock access to the “public good” of new space via non-cell autonomous action of environmental engineering, mediated by the action of secreted ECM-remodeling enzymes. We found that non-cell autonomous production of common good had a profound impact on both phenotypic and clonal evolution, leading to higher heterogeneity and the emergence of patterns of phenotypic succession, while at the same time producing clonal architecture that might be consistent with a hypothesis of evolution neutrality.

Consistent with published reports of convergent evolution in cancers, where mutations with similar functional consequences are observed in distinct subclonal sub-populations (Ardaševa et al., 2020), we observed phenotypic convergence in our simulations, as different lineages within computational tumor were capable of finding optimal phenotypic solutions (Gerlinger et al., 2012). Although consideration of non-cell autonomous environmental engineering increased phenotypic diversity, it did not qualitatively impact convergence.

The addition of public good-enabled populations of tumor cells led to higher levels of clonal heterogeneity. To a large extent, this effect likely reflected the ability of these tumors to reach a much higher population size, as it was most pronounced within smaller subdomain size. However, public good engineering action also enabled a more effective maintenance of clonal ITH (notice the lack of a dip after reaching the peak in green versus black lines in Figure 2A) – consistent with the report of maintenance of clonal heterogeneity in a mouse model of non-cell autonomous driver of tumor growth (Marusyk et al., 2014).

The inclusion of producers of public good leads to more complex ecological dynamics, involving local and global phenotypic successions. The ability to degrade ECM barriers can, in some cases, lead to strong selection for engineering phenotypes, as they were capable of getting priority access to the newly available space – making the engineering strategy more effective at the growing edge of a tumor. However, as proliferative phenotypes were more successful in utilizing space, they eventually outcompeted engineers, leading to robustly reproducible patterns of phenotypic succession, which parallels ecological succession.

Even though the dynamic behavior in our simulations was designed to be governed by strong selection for competing phenotypic strategies, it produced patterns of clonal architecture that could be interpreted as indicative of neutral evolution patterns within current approaches toward inferences of evolution mode from genomic sequencing data. Most likely, this result reflects a lack of mutational record of history of phenotypic switching. Surprisingly, inferences of neutrality based on goodness of fit between the inverse of VAF and cumulative numbers of mutations showed higher consistency of the outcomes of our selection-driven simulations compared to simulations with explicitly neutral dynamics. Although increased resolution of mutation sampling led to superior fit in neutral simulations, the effect was only observed at resolution levels that exceed those employed in typical mutational analyses studies, with the exception of analyses of single crypts in colorectal cancers. Although our results do not disprove recent reports of neutral evolution in many cancers (Williams et al, 2016, 2018), they do warrant a more cautious interpretation of evolutionary inferences based on mutational analyses.

Limitations of the study

Although mathematical and computational modeling enables quantitative analysis of evolutionary dynamics in space in time, modeling deals with abstraction and simplification of biological reality. Moreover, the use of computational modeling depends on the accuracy of model assumptions and parameters, which are specified explicitly. While tumor growth requires gaining access to new space, this access is mediated not only by the action of tumor cells but also by complex networks of microenvironmental interactions involving non-tumor cells, such as fibroblasts, macrophages, etc (Janiszewska et al., 2019; Marusyk et al., 2014). A concrete understanding of the underlying biology at this level, permitting adequate capture in a mathematical model, is still missing. Even if such knowledge were available, quantitatively specifying all of relevant interactions would not be feasible. While primary tumors might be simpler than natural ecosystems, their ecological complexity is certainly higher than what can be captured in our models (Brown, 2016). Further, two-dimensional lattices considered in our work do not faithfully recapitulate three-dimensional tissue structures within primary or experimental tumors. On the other hand, computational modeling enables the examination of outcomes when the key ecological and evolutionary drivers of population dynamics are considered – something which is not accessible not only to inferences from analyses of clinical samples but also to experimental studies. Therefore, while quantitative results of our simulations might not be directly applicable to understanding of tumor biology, qualitative finding of the impact of ecological engineering on tumor heterogeneity and phenotypic succession is likely to be relevant to real tumors. Non-cell autonomous actions, including public goods, appear to be common in tumor ecosystems, including access to vasculature nutrients, immune invasion, etc (Tabassum and Polyak, 2015). Although our simulations only considered non-cell autonomous effects in terms of accessing spatial resources, results from our simulations are likely to be generalizable to other scenarios involving other public goods. We posit that consideration of ecological aspects of cancer evolution, the presence of non-cell autonomous effects, as well as integration of mathematical modeling with experimental studies, and analyses of clinical samples will be essential to better understand clonal evolution – which is a prerequisite to the development of successful strategies to prevent and eradicate advanced cancers.

Resource availability

Lead contact

Further information and requests for resources and data should be directed to and will be fulfilled by the Lead Contact, Jack Edwards (jackedwards1@usf.edu).

Materials availability

This study did not utilize any physical materials nor does it contain any biological data.

Data and code availability

The model is publicly available on GitHub at CE_ABM: https://github.com/jackedwards1/CE_ABM.

The accession number for the data used to produce the figures reported in this paper is available on Mendeley: https://dx.doi.org/10.17632/5wpcrnkvc7.1.

Methods

All methods can be found in the accompanying Transparent Methods supplemental file.