C ancer evolution in a changing microenvironment

Cancer development can be viewed as an evolutionary and ecological process, in which the tumour microenvironment (TME) is likely to play a critical role. Unfortunately, the TME is largely ignored or considered static in most cancer evolution models, and different cancers are often studied in isolation. A general theory of adaptive cancer evolution is lacking. Here we establish a phenotypic and genetic model of cancer evolution in three-dimensional (3D) space with a changing TME. With individual-based simulations, we show how cancer cells adapt to diverse changing TME conditions and fitness landscapes. Compared with static TMEs, changing TMEs can generate complex 3D dynamics of spatio-temporal heterogeneity involving variable subclonal fitness and mixing, driver mutations with different fitness effects and phylogenetic patterns. Our 3D simulations with changing TMEs capture some of the key morphological characteristics of cancer, including spatio-temporal ball-like and irregular clonal/subclonal population structures. A cycling TME, in particular, is capable of generating more driver mutations and promoting cancer adaptation. We predict that the TME is a major limiting factor of adaptive cancer evolution. Finally, our model can be used to simulate anti-cancer treatment strategies and show how they can be subverted by different resistance mechanisms. Our study provides an evolutionary and ecological framework for understanding cancer development and treatment, and provides novel insights into the processes of adaptive cancer evolution and precision cancer medicine.


Introduction
Understanding the mechanisms facilitating ecological adaptation is a fundamental question in evolutionary biology, including cancer development 1 . In 1976, Nowell outlined carcinogenesis as an evolutionary process using an illustrated phylogenetic model in which clonal and stepwise accumulation of advantageous mutations under selection from the tumour microenvironment (TME) was envisaged 2 . Our understanding of cancer development as an evolutionary process has increased greatly in the last decade due to advances in next generation sequencing and related analytical methods 3 .
Many ideas and concepts originating in evolutionary biology and molecular evolution have become widespread in cancer research including gradual evolution with stepwise accumulation of driver mutations 4 , punctuated evolution with "catastrophic" genomic changes 5 and neutral cancer evolution 6 . In some cases, anti-cancer treatment has been developed based on evolutionary principles 7 . Many interesting patterns have been observed based on single time point data from cancers at their time of surgical resection, such as spatial and temporal heterogeneity of subclonal mutations, subclonal mixing 8,9 , and ball-and non-ball-like clonal and/or subclonal structures 10,11 .
If the environment plays an important role in determining which somatic mutations drive tumour evolution under selection, the cancer cell and its TME should jointly determine the tumour evolutionary trajectory. Whilst the cancer cell may determine its own microenvironment to some extent, the source of selection extends to include many non-neoplastic cells, including various stromal components such as fibroblasts, immune cells, the extracellular matrix (ECM) and the pre-metastatic niche 12 . The TME is thought to be a dynamic regulator of tumour progression and metastasis 11,[13][14][15][16][17] . Finally, the origin of cancer cells has been extended to both stem cells and non-stem cells 18,19 .
Although the seminal ideas proposed by Nowell and others -that tumorigenesis proceeds by accumulating stepwise advantageous mutations under selection -have guided cancer research for decades, there is a lack of a general theory of adaptive cancer evolution. Some fundamental questions remain unanswered, including the role of the (changing) TME in determining the fitness effect of mutations and the evolutionary trajectories of cancers in three-dimensional (3D) space. Does cancer evolution proceed by mutations with small fitness effects, large effects or a combination of both? How and why does the fitness effect of mutations differ across cancer types? Finally, what spatiotemporal patterns of cancer development can be observed under different tempos and modes of adaptive cancer evolution?
In previous cancer evolutionary modelling, in most cases the TME has been neglected and the fitness effect of mutations must therefore be pre-specified 6,10 . When the TME has specifically been considered as a source of selection, the underlying genetic basis [20][21][22] , such as the variable fitness effects of adaptive mutations, have not been explicitly considered. Moreover, the important question of how cancer cells adapt in a changing TME has not generally been considered. To address these issues, we have analysed carcinogenesis as an adaptive evolutionary process, in which the cancer cell and its changing TME jointly determine the cancer evolutionary trajectory in 3D space, with the aim of providing a general theory of adaptive cancer evolution (see Methods for details).
By simulating cancer growth in the 3D space of a changing TME in different fitness landscapes at single cell and population level, we identify various evolutionary patterns of 3D adaptive cancer evolution under diverse genetic, phenotypic, population genetic and changing TME conditions. We show that far more complex evolutionary patterns can emerge compared with models that assume a static TME. Important differences include the spatio-temporal heterogeneity in clonal/subclonal fitness, population mixing and structures, and the combination of driver mutations with various selective advantages 23 .

Description of the model
Our model uses formal adaptation theory and simple 3D growth dynamics, based on R.A.
Fisher's geometric framework. This was originally described as the phenotypic geometric model in his seminal 1930 work "The genetical theory of natural selection", which provided much of our understanding of the nature of adaptation inherent in Darwin's theory of evolution by natural selection. The adaptive process Fisher described is analogous to adjusting the knobs of a microscope and each trait contributes to fitness independently (universal pleiotropy) 24 . We assume that a changing TME in 3D leads to continual selection acting on the proliferation potential (the fitness effect) of mutations that influence the phenotypic traits of cancer cells. This is modelled by cancer fitness landscapes, with a single changing phenotypic optimum ( v 1 ) determined by the TME (see Methods and Supplementary Figure S1).
In our model, each cancer cell is, for simplicity, diploid. We consider n traits and L loci.
The number of traits is related to the notional phenotypic complexity of the cancer cell.
It also determines the dimensionality of the fitness landscape, with a higher number of traits indicating a more "complex" cancer cell (for example, including more of the hallmarks of cancer 14 , with traits such as aerobic glycolysis and immune evasion). In a very simple situation, the fitness of a cancer cell with two adaptive traits can be described by the fitness function in a three-dimensional Cartesian co-ordinate system, where the height along the fitness surface corresponds to fitness and the other two coordinates (Supplementary Figure S2) correspond to phenotypic values of each trait. The initial optimum phenotypic value defining the maximum fitness (the "peak" of the fitness landscape, see Supplementary Figure S2) is at the origin of the Cartesian coordinate system, and the optimum can move, directionally or otherwise (e.g., randomly or cyclically, see Methods), as the TME changes.
For a cancer cell, a phenotype, z , of n traits is defined by a column vector, z = z 1 ,...,z n ( ) T (equation (1)). We assume that the optimum, z opt = z 1 opt ,...,z n opt ( ) T (equation (2)) moves away from the origin following different dynamics (see Methods). So for a particular point representing a cancer cell's phenotype and fitness in the fitness landscape the Euclidean distance d between phenotype z and its optimum z opt (equation (3)) and the fitness function w(d)(equation (4)) determines the cancer cell's fitness (if epistasis is not considered). We assume each mutation affects all n traits according to the assumption in Fisher's geometric model. Each trait can then be changed with a certain probability by a random mutation, defined by a column vector, r = r 1 ,...,r n ( ) T (equation (5)) that represents the phenotypic effect of the mutation at any of the 2L loci (see Methods for details). The mutations in each generation additively determine the new ndimensional phenotype ′ z , which can then be changed again by further random mutations.
Both the moving phenotypic optimum of the TME and random mutations can change a cancer cell's fitness. For example, by equation (4), when the optimum moves away, the relative distance, d , will increase to d ′ , which leads to a smaller fitness value w( ′ d ).
Similarly, a random mutation can also decrease or increase the phenotypic Euclidean distance to ′ d and lead to a new fitness w( ′ d ), which may be higher or lower depending on the phenotypic effect of the random mutation. The phenotypic effect (size) of a mutation is sampled from a multivariate normal distribution (equation 11). A mutation with large fitness effect can move the cancer cell a long phenotypic distance relative to its phenotypic optimum in the fitness landscape. All these changes in fitness due to either mutational or TME effects naturally lead to different levels of selection and adaptation.
At each generation, a cancer cell is selected according to its fitness value and dies with probability 1− w z,t ( ) , which depends on the phenotypic effects of the mutations in that cell and position of the TME optimum at time t. The surviving cancer cells reproduce asexually. New daughter cells randomly occupy available space in 3D according to prespecified limits on tumour size. The tumour may go extinct, persist for a long time with varying size, or grow continuously until it reaches its maximum allowed space or size, after which it continues to be under viability selection without expanding. A cancer cell acquires higher fitness when mutation(s) move its phenotype to be closer to the We use computer simulations based on the above principles to assess cancer growth in 3D and look for patterns of adaptive cancer evolution. The initial fitness ( w 0 = w(z 0 )) of the cancer cell is conferred by a starting phenotype, z 0 , which is pre-specified and could be determined by a new driver mutation or conceivably by phenotypic plasticity, which may, for example, arise from disturbances such as inflammation 13 . Mutation rates and the number of loci are also determined before each simulation: the mutation rate could, for example, be increased owing to mutations in key pathways maintaining DNA replication fidelity, or extrinsic mutagens; and the number of loci could also be increased or decreased, for example owing to changes in somatic copy number variation.
In the following sections, we first show how cancers adapt under the classic model of tumorigenesis, where the tumour is initiated from a single stem cell with high fitness, and we identify the general patterns of 3D spatio-temporal and adaptive cancer evolution in a dynamic TME with different rates of change. We then simulate cancer evolution following more general assumptions of tumour cell origin, such as different initial population sizes, phenotypes, fitness, and shapes of the fitness landscapes and other dynamics of TME change. Finally we show that our model can be used to demonstrate the effects of anti-cancer treatment strategies, such as radiotherapy, chemotherapy and immunotherapy, where the treatments can cause sudden or gradual changes in the TME optimum, thus changing the selection dynamics; we also show how different resistance mechanisms can subvert an apparently effective treatment.
Comparing adaptive cancer evolution in three dimensions under various changing and static TMEs For simplicity we first assess the classic model of tumorigenesis in a directionally changing TME. A single cell with two traits ( n = 2 ) starts asexual reproduction from the centre of the 3D tumour space with a phenotype at the origin with initial fitness w 0 = 1 , . CC-BY-NC-ND 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/377226 doi: bioRxiv preprint first posted online Jul. 25, 2018; which could represent a stem cell with high fitness conferred by a driver mutation that starts to form a tumour in its local TME. The cancer starts to grow with cell birth probability, w z,t ( ) , and death probability, 1− w z,t ( ) . The cancer may go extinct if it has low average fitness, while higher average fitness can keep the cancer persistent for long time.
In a static TME ( v 1 = 0 , the optimum does not move), the cancer can persist as long as an Second, we find that the speed of TME change plays a critical role in cancer progression and adaptation. Particularly, if the TME optimum moves too fast and the cancer cells do not acquire enough beneficial mutations to increase fitness, the population goes extinct . CC-BY-NC-ND 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/377226 doi: bioRxiv preprint first posted online Jul. 25, 2018; Third, we find two general clonal or sub-clonal growth patterns in 3D: (i) one or more balls of cells (Figure 1a-i); or (ii) irregular morphology (Figure 1j-l). Intriguingly, these patterns generally resemble those observed clinically 10,11 . It is more likely that ball-like clones or sub-clones are found in a fast-changing TME ( v 1 = 0.05 and v 1 = 5×10 −3 , Figure   1a-f, Supplementary Movies S3 and S4) than a slow-changing TME ( v 1 = 5×10 −4 and v 1 = 5×10 −5 , Figure 1g-l, Supplementary Movies S5 and S6). Moreover, the fast-changing TMEs also lead to many smaller sub-clones, which turn over fast, with fluctuating numbers, spatial proximity and fitness (Figure 1a We then sought to understand the fitness effects of mutations fixed during cancer development under a changing TME. First, in a static TME, although the cancer can persist for an individual's natural life, we find that there is no fixation of any new driver mutations (Supplementary Figure S3 and Movies S1). Interestingly, this is a more general finding that might apply to the "Big Bang" model of colorectal tumor growth, which is characterised by a lack of selective sweeps and high intratumoural heterogeneity 6 . In a changing TME the mean population fitness and adaptive mutations are determined by the speed of TME change ( Figure 3). As expected, the mean population fitness decreases faster in a fast-changing TME, whereas a slowchanging/static TME leads to slower fitness decay and thus longer survival time of the tumour ( Figure 3a). Lower selection intensities (a "flatter" fitness landscape; Supplementary Figure S2 Figures S4 and S5).
Second, in a fast-changing TME, the driver mutations that spread to fixation must have in general larger fitness effects (large selection coefficients, Figure 3b). This can be . CC-BY-NC-ND 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/377226 doi: bioRxiv preprint first posted online Jul. 25, 2018; explained from equations (12) and (13), as for an individual cancer cell, a faster changing optimum will produce relatively longer distances from its location to the optimum in the fitness landscape, thus reducing fitness and meaning that the cell is less likely to pass viability selection. In order to avoid extinction, cancer cells must continuously acquire mutations with larger fitness effects to catch up with the fast moving optimum, which becomes impossible when the TME changes too fast. These cancers can only survive for relatively short periods of time and harbour a few large adaptive mutations (Figure 3b Intriguingly, a moderately changing TME promotes cancer evolution with the highest  Figure S4). This type of evolution by advantageous "nearly neutral" mutations with small fitness effect, as Tomoko Ohta pointed out, could be interpreted as adaptive evolution in a changing environment 26 .
Finally, when the number of traits of cancer cells under selection increases from 1 to 8 (increase of cancer cell phenotypic complexity, n , see Methods), we find that there is a cost of complexity associated with cancer adaptive evolution. In particular, the mean fitness of the population decreases and the mean selection coefficient increases during adaptation when the cancer cell has an increased number of traits and changing TME (Figure 3d-e). This observation could be understood in terms of the hallmark traits of cancer cells 14 . Suppose that a tumour cell population "devotes" all its genetic variation to the adaptation of one or two key traits instead of eight traits, adaptation may become more efficient because there is no requirement for larger driver mutations. As a result . CC-BY-NC-ND 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/377226 doi: bioRxiv preprint first posted online Jul. 25, 2018; the cancer with lower phenotypic complexity may persist for longer period of time and it is more likely to establish a clinically significant cancer.
3D adaptive cancer evolution under diverse cancer initiation conditions and a changing TME As stated above, a single cell with initial high fitness ( w 0 = 1 ) due to driver mutations is almost guaranteed to initiate cancer growth under a directionally moving TME optimum.
However, a cancer cell with lower fitness ( w 0 < 1) may also initiate tumour growth under different TME changing directions and rates. To gain further insights, we now assume more general assumptions for cancer initiation and TME changing dynamics, such that the neoplastic growth is initiated by somatic cells with different initial phenotype/fitness and population sizes 19 . This can also serve as a model of metastasis in that migrating cells with different phenotypes and fitness must adapt at a distant site with a novel TME. In this general model, 3D cancer growth can in principle be initiated from any number of cells with any fitness in a TME that can take different types of changing dynamics.
We consider that the population starts from a phenotype away from the optimum, and thus has low initial fitness (equation (9), w 0 < 1). In order to avoid extinction (Supplementary Movie S7), a cancer must be initiated with a large population size when its initial fitness is very small (e.g., w 0 = 0.1 ) and initial driver mutations must occur with relatively large fitness effect (Figure 4a-d, Supplementary Movies S8-S11), which is true for both static and changing TMEs. In a static TME ( v 1 = 0 ) when the initial cancer cell is away from the optimum (e.g., w 0 = 0.6 ) our simulation shows that only two driver mutations are recorded in subsequent adaptation (data not shown, similar to Figure 4e). This type of mutation should be considered to be classical drivers. This means that low fitness cells (a long phenotypic distance from the optimum) with a large initial population size have a higher chance of generating appropriate driver mutations to survive initial strong selection due to either sudden or rapid TME change. Higher fitness cells require few cells to initiate the neoplastic growth, which may be conferred by . CC-BY-NC-ND 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/377226 doi: bioRxiv preprint first posted online Jul. 25, 2018; either phenotypic plasticity or classical driver mutations already present in the initial transformed cells (movies not shown for w 0 = 0.9 as they evolve similarly to w 0 = 1 , see When the population survives its initial selection with driver mutations, its evolutionary trajectory still depends on the rate of TME change (Figure 4, Supplementary Movies S8-S15). Again our simulations show two similar patterns as above. First, a moderately changing TME promotes cancer adaptation by fixing more driver mutations (Figure 4b, f and j, Supplementary Movies S9, S13). Second, if the population evolves under a slowly changing TME, we recover subsequent driver mutations that have small fitness effects (Figure 4a, e and i, Supplementary Movies S7, S15). These mutations could be termed "mini-drivers", as we proposed previously 23 . Interestingly, this phenomenon resembles the pattern of Fisher's micromutationism, in which adaptation proceeds initially by mutations of larger fitness effect and subsequently smaller fitness effect mutations 27 predominate when the phenotype approaches the optimum (assuming that the TME optimum moves quite slowly, Figure 4).
In cancer evolution, most driver mutations have pleiotropic effects, affecting multiple pathways and several cancer traits. Although we assumed that each mutation affects all traits (universal pleiotropy), its fitness effect on each trait may be different. In other words, the driver mutations may have different fitness effects on different traits and the selection is therefore correlated (equation (9)), which can be illustrated by the shapes of the fitness landscapes with different levels of correlation (see Supplementary Figure S2).
Similarly, a mutation (genotype) may have different phenotypic effects on each trait. So the mutational effects on traits can be independent or correlated (see mutational distribution, equation (11)). Indeed, we find that when the TME changes slowly, the fitness effects of driver mutations along traits indeed correlate (Supplementary Figure   S6). Moreover, we find that when the TME changes fast the phenotypic effects of driver mutations along traits have similar levels of correlation as correlations of mutational effects (data not shown). Although in clinical settings the effects of selectional and mutational correlations on cancer progression are unclear, an example may be multiple . CC-BY-NC-ND 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
Finally, to understand how other TME changing dynamics may affect cancer adaptation we analysed another three types of changing TME. First, when the TME changes randomly, the resulting increase in the TME variance at different time points always acts against cancer adaptation (equations (14)-(16)), leading to reduced mean population fitness and a requirement for adaptive mutations with higher mean fitness effect (Supplementary Figure S7). Second, when we add a random component into a directionally changing TME (equations (17)-(19)), the increased variance caused by the random component also acts against cancer adaptation, producing similar results to a purely randomly changing TME (Supplementary Figure S8). Third, when the TME changes cyclically (equations (20)-(22)), with increased amplitude, the cycling TME optimum also acts against cancer adaptation (Supplementary Figure S9). However, interestingly, although the mean population fitness decreases and mean selection coefficient of adaptive mutations increase when the amplitude increases, there are more adaptive mutations recorded than under any other TME changing dynamics (full data not shown, amplitude A = 4 , see Supplementary Figure S9 and Supplementary Movies S16-S17). Moreover, there are also more complex spatio-temporal patterns of sub-clonal fitness and mixing, in which birth and death of large and small sub-clones with diverse fitness values occur frequently through time and space (Supplementary Movies S16-S17). This indicates that a cycling TME at intermediate level may be particularly capable of promoting cancer adaptation by periodically fixing more adaptive mutations than TMEs that change directionally and/or randomly.
We then reconstructed cancer phylogenies in the three types of simulations. Intriguingly, we found that in all cases the shape and the temporal signals in the phylogenies (e.g., the branch length and overall shape) were characteristic of its particular TME changing dynamics, such as static, directional, random and cyclic TMEs with different optimum changing speed (illustrated and explained in Supplementary Figure S10). These patterns from cancer phylogenies could be particularly useful in inferring the underlying TME selection dynamics and cancer evolutionary history.
. CC-BY-NC-ND 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/377226 doi: bioRxiv preprint first posted online Jul. 25, 2018; Cancer adaptation owing to anti-cancer therapies Anti-cancer therapies of all sorts, but principally those that are non-surgical, introduce an artificial TME that can radically change the tumour's selective landscape, at least transiently and possibly permanently. In the case of targeted therapies, for example, cells with a beneficial mutation can rapidly become at a selective disadvantage. As a proof of concept, we find that our model can be used to understand anti-cancer therapies that use different dosing strategies, by assuming that these modify the TME selective optimum and/or the shape of the fitness landscape (see equations (25)-(28), Supplementary Figure S2). Initially, we have examined treatments that cause a sudden change in the TME optimum, such as genotoxic therapies (radiotherapy or chemotherapy) 28 and immunotherapy 29 (e.g. immune checkpoint inhibitors), which should put strong selective pressure on all cancer cells.
For simplicity, we assume that the therapy causes the TME optimum to move suddenly along the axis of a single trait. We first assume that the cancer has evolved to a constant TME optimum z 0 opt = 0 and allow it 100 generations to accumulate genetic variation and reach maximum tumour size ( N = 1×10 7 ). We then assess four different TME optima to represent different treatment strategies ( z 1 opt = 5,6,7,8 , see equations (25)-(28), Supplementary Movies S18-S21) that treat the cancer for about 33 months or more (e.g., 1000 generations), reducing mean population fitness below w = 0.1 . We illustrate how the optimum of the fitness landscape changes from z 0 opt to z 1 opt under each treatment (Supplementary Figure S11).
As shown in Figure 5, all treatments reduce the fitness of all cancers below w = 0.1 (Figure 5a-o) and at z 1 opt = 8 the cancer is successfully cured (the population is extinct, Figure 5p), but due to mutation and smaller TME change (e.g., because of smaller dose or effectiveness of the delivery) the cancers survive at z 1 opt = 5,6,7 and quickly relapse in less than two months (Figure 5p). In general, when the dose is not sufficient for maximum killing, smaller effective doses lead to selection for mutations with smaller . CC-BY-NC-ND 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/377226 doi: bioRxiv preprint first posted online Jul. 25, 2018; fitness effects, while higher doses lead to the opposite. Phylogenetic trees show that treatment selection often leads to expansion of sub-clones and extinction of the majority of the tumour mass (Figure 5e, j and o). Weak treatment z 1 opt = 5 leads to a large number of initial sub-clones with resistant beneficial mutations, but eventually the tumour mass comes to be dominated by a single sub-clone (Figure 5e) carrying the mutation with the largest fitness effect ( s = 5.19 ). Stronger, but non-lethal, treatment ( z 1 opt = 7 ) leads to an early sub-clonal expansion carrying the mutation with a very large selection coefficient ( Figure 5o).
In terms of sub-clonal diversity, we first observe a reduction soon after treatment and then an increase of diversity, but after the sub-clonal fitness recovers to the same/similar levels as pre-treatment or with minor loss, the diversity reduces again until the end of the treatment (data not shown). These behaviours are consistent with a mean fitness decrease due to the treatment and selection for resistant sub-clones ( Figure 5p). Interestingly, our simulations show that both pre-existing and de novo resistant mutations can be observed. Moreover, when the treatments lead to higher selection intensity (a "narrower" fitness peak, see Supplementary Figure S2), all treatments lead to immediate cure (population extinction).
To understand potential therapeutic resistance, we examined three mechanisms that We now combine our findings from our general model of cancer evolution with those of the evolutionary responses to anti-cancer therapy. Treatment strategies that lead to a moving TME optimum (in any direction) may be effective and help to reduce toxicity by using of a lower dose, as there is no requirement of initial maximum dosing to induce a high sudden optimum change for maximum killing ( in directional dosing the treatment must be stopped when the allowed dose is reached or continue if the tolerance is increased. Since moving TME-based anti-cancer strategies are highly desirable, designing clinically effective dosing and delivery strategies is important, although currently challenging.

Discussion
In this study, we provide an evolutionary and ecological framework for understanding cancer development and anti-cancer treatment strategies. This framework captures the complex 3D spatio-temporal dynamics of intra-tumour heterogeneity in sub-clonal fitness and structure, and the tempo and mode of adaptive cancer evolution in a changing TME. Cancer evolution appears to be far more complex in a changing TME (e.g., . CC-BY-NC-ND 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/377226 doi: bioRxiv preprint first posted online Jul. 25, 2018; Supplementary Movies S1-S6, S16-S17), particularly as regards the early stages, including rapid sub-clonal shifts and "spontaneous" regression or tumour death, which suggests there are many unobservable failed evolutionary trajectories of cancer progression in human patients unless caught prior to extinction 17 . Our model is thus biologically interpretable and intuitive, as parameters of birth, death and selection are naturally incorporated into the properties of the TME and the underlying phenotypic effect of mutations.
Currently the actual fitness effects of most mutations in cancer genomes are unknown, although estimates suggest that selection is often weak, even for major driver mutations 30 . Our model predicts that there are many mutations with small fitness effects due to weak selection caused by a slowly or moderately changing TME. These mutations may be pervasive in current cancer genomic data but are difficult to identify due to their small effects. Our study gives a theoretical basis to the observed weak selection, and may help to explain why cancers are rarer than would be predicted given the number of cells in the body and the potential role of the TME in restraining cancer progression 17 .
Many cancer driver mutations might be classified as "nearly neutral" 31 as described by Tomoko Ohta 26 . However, as Tomoko Ohta pointed out, due to a changing environment, it is appropriate to classify fixed advantageous mutations as positively selected, even when their fitness effect is quite small. So it is appropriate to term these mutations as "mini drivers" 23 .
Adaptive cancer evolution in individuals may have several evolutionary tempos and modes mixed across cancer types in humans, which can be explained by the heterogeneity of the TME and its varying optimum. On one hand, cancer cells may frequently go to extinction due to strong stabilizing selection from the normal TME. On the other hand, an extremely slow-changing or constant TME may lead, in effect, to limited cancer cell evolution in which neutral/nearly neutral mutations accumulate, as might be the case for cancers that apparently carry no or few classical driver mutations.
Of course, the TME may determine the phenotypic effect of a mutation as well as determine the selective landscape, so a changing TME might alter both a phenotype and its selective advantage. Furthermore, cancer growth and progression (e.g. acquisition of new mutations) will also affect the TME; and the cancer itself may create or modify a . CC-BY-NC-ND 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/377226 doi: bioRxiv preprint first posted online Jul. 25, 2018; TME optimum leading to non-cell-autonomous cancer evolution. Therefore, quantitative understanding of these basic components and TME changing dynamics in individual patients is crucial for developing future cancer medicine. For instance, we show that a cycling TME may be particularly capable of promoting cancer adaptation (see Supplementary Movies S16-S17 for its unusual long-term sub-clonal dynamics and a large number of recorded adaptive steps), while stochasticity and fast changes in the TME may act against cancer adaptation. TME dynamics need to be considered alongside other factors in studies that aim to provide a complete picture of the underlying cancer evolutionary dynamics.
In this study we show that in a moderately/slowly changing/static TME, lifetime seeding is possible to distant sites (e.g., see Figure 3a and 3b, v 1 = 5×10 −4 and v 1 = 5×10 −5 ), but successful colonization at distant metastatic sites is stochastic and uncertain, because it depends on the initial phenotype and population size of the circulating cancer cells and their interactions with the distant metastatic TMEs (or pre-metastatic niches) assuming that metastatic potential depends on the long survival of the primary cancer due to a slowly changing TME (Figure 4). Phenotypic plasticity may play a role in initial adaptation, which confers higher fitness of initial cells with a small migrating population (Figure 4i-l). This is consistent with the hypothesis that the pre-metastatic niche may bring the phenotype of the migrating cancer cells closer to the optimum of a potential metastatic TME for higher fitness 12,15,32 . There is empirical evidence that cancer cells migrate in small numbers. Metastasis does not necessarily require that circulating cancer cells have major driver mutations to initiate cancer adaptation 33 as shown in Figure 4i. Although in a static TME our observation does support the main characteristics of the "Big Bang" model, we show that the main evolutionary force at play in this situation is likely to be weak purifying selection 34 . Our simulations may also capture some of the unobservable aspects of cancer development, which may happen before the "Big Bang".
Current cancer data are mostly derived from relatively late stages of cancer development, where the cancer has accumulated large heterogeneity under weak dynamic selection possibly due to a changing TME. Although the initial classic driver . CC-BY-NC-ND 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/377226 doi: bioRxiv preprint first posted online Jul. 25, 2018; mutations with large fitness effect can be easily detected using next generation sequencing, the fitness effect of such mutations under a changing TME remains unclear in clinical settings, which leads to scenarios where no obvious drivers can be detected with current approach when the fitness effect of these driver mutations is small. This poses significant challenges for cancer treatment. Our results support the use of proposed therapeutic strategies that can target the TME, causing to change, and drive cancer cell populations to extinction 17,35 . Importantly, the TME dynamics required to kill cancer cells may be different for each individual patient, which will further require precise and individualised modelling of therapeutic dosage and delivery strategies.
Although a similar approach, termed "adaptive therapy", has been proposed previously 36 , where treatments have been adjusted by evolutionary principles and there is clear evidence of improved patient outcome by cyclic dosing 37 , the true evolutionary dynamics of cancer progression and treatment can only be understood by combined measures of cancer genetic and phenotypic changes and the corresponding TME dynamics, as we demonstrate in this study.
In conclusion, our 3D model provides a natural evolutionary and ecological framework to understand adaptive cancer evolution. We show that the extension of Fisher's classical model using fitness landscapes with a changing TME is sufficient to produce many of the observed complex cancer evolutionary and 3D pathology patterns. In the future, in vivo and/or in vitro experiments, combined with such modelling approaches and fine-scale methods, such as single-cell sequencing and phylogenetics, could be used to elucidate cancer fitness landscapes with dynamic TMEs, further testing the (un)predictability of cancer evolution, with implications for cancer precision medicine 38 .

Model
In evolution, organisms can evolve in a Darwinian or non-Darwinian way or a combination of both depending on how natural selection acts on the phenotypic traits and their plasticity 39 . The fitness advantage to organisms conferred by these traits due to environmental selection can lead to adaptation [40][41][42] . However, adaptation may not . CC-BY-NC-ND 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
For instance, a plausible path of adaption, which involves both plasticity and selection, starts with the initial plastic adaptation (initial fitness advantage) and subsequent stabilizing selection and/or positive selection due to environmental changes 15,39,46 .
These evolutionary theories help to define a general mode of cancer initiation and subsequent adaptation. In this study, we follow Fisher's geometric framework 27,[47][48][49]  The traits involved can be any, but could for convenience be those highlighted by Hanahan and Weinberg 14 : sustaining proliferative signaling, evading growth suppressors, resisting cell death, enabling replicative immortality, angiogenesis, activating invasion and metastasis, reprogramming of energy metabolism and evading immune destruction 14 . The number of traits or the dimension of the phenotypic vector, n , represents the "complexity" of a cancer cell. As we assume that the TME (the ecology) is the primary source of selection 17,51-56 , there is a corresponding optimum phenotype z opt , which is defined by a column vector of n values The Euclidean distance, d , is defined as . CC-BY-NC-ND 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/377226 doi: bioRxiv preprint first posted online Jul. 25, 2018; The phenotypic fitness function for a cancer cell in a microenvironment is defined as 49 Here a is the selection intensity for all traits ( a > 0 ) 49 . Therefore, the fitness of an individual cancer cell depends on its phenotype's Euclidean distance to the optimum.
Equation (4) suggests that the closer a cancer cell's phenotype is to the optimum, the fitter it becomes. The change of phenotypic traits of a cancer cell due to random mutations can be defined as an n-dimensional random number: the size of a mutation (the effect of a mutation on phenotypic traits) is therefore defined so the combined phenotypic trait of a cancer cell ′ z with mutation r , relative to its wild type z is defined as . CC-BY-NC-ND 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/377226 doi: bioRxiv preprint first posted online Jul. 25, 2018; the selection coefficient of a mutation that changes the fitness of the cancer cell is therefore defined as when s > 0 the mutation is beneficial, it moves the cancer cell closer to the optimum.
When the cancer cell is at optimum we say it has its maximum fitness w z ( ) = 1 , whereas when the cancer cell moves away from the optimum its fitness decreases and it could reach its lowest fitness w z ( ) = 0 . Such mutations are deleterious and lead to negative selection coefficients with s < 0 . When mutations do not change fitness, they are defined as neutral mutations with s = 0 .
We can now define a general form of the fitness function from equation (4) as shown before 24,47,48,57 .
where S is a real n× n positive definite and symmetrical matrix, and T denotes transposition. The matrix S describes the shape of the fitness landscape, namely, the selection intensity. Matrix S −1 is the inverse of S . As introduced above, the overlapping of pathways responsible for cancer cell traits indicates pleiotropic effect of mutations . CC-BY-NC-ND 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/377226 doi: bioRxiv preprint first posted online Jul. 25, 2018; contributing to cancer adaptive evolution 14 . If the selection intensity is the same along all n traits then we have an isotropic fitness landscape (universal pleiotropy). We set S = σ 2 I ( σ 2 > 0 , I is an identity matrix). Selection may also vary along different traits (selection is correlated), which means mutational effects contribute to fitness differently for different traits. If we set selection intensity to vary along n traits, then S has nonzero off-diagonal entries. We illustrate different shapes of the fitness landscape (independent and correlated selection, Supplementary Figure S2). We can use to measure the average width of the fitness landscape, where σ 2 is defined as the geometric mean of the eigenvalues of S . The selection coefficient of a mutation in equation (8) can now be defined as

Distribution of mutation effect size
We assume the phenotypic effect of a mutation with size follows a multivariate normal distribution with mean 0 and n× n variance-covariance matrix M 48,58: as above M is also symmetrical and positive-definite. We have M = m 2 I , where m 2 is the variance of mutational effects and I is an identity matrix. If mutational effects are correlated we have covariance of mutational effects as the off diagonal elements in M .
r . CC-BY-NC-ND 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/377226 doi: bioRxiv preprint first posted online Jul. 25, 2018; The average variance of mutational effects, m 2 , is given by the geometric mean of the eigenvalues of M . We have m 2 = det M ( ) n .

Cancer adaptation in a changing microenvironment
In all, we consider six scenarios for the properties of the tumour microenvironment.
Different initial population or individual cell properties (e.g., population size and fitness) during computer simulations are specified in the results section. We assume that the TME optimum of the first trait of the n traits changes 48,50 .

Directionally changing optimum
We assume a directionally moving optimum with vt where the vector v = v 1 ,...,v n ( ) ′ is the TME change speed. The optimum z opt t ( ) is time (generation) dependent and only the first trait optimum changes.

Randomly changing optimum
Here we assume the TME optimum of the first trait, z 1 opt , changes randomly at each generation, t , following a Normal distribution with mean 0 and standard deviation, δ z opt t ( ) = z 1 opt ,...,z n opt

Directionally changing optimum, with a random component
We assume a directionally moving optimum with vt (see equation (12)) and a random component e . CC-BY-NC-ND 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

Cyclically changing optimum
In this scenario, we assume the TME optimum changes cyclically or periodically at each generation, t , following a general periodic function 50,58,59 z opt t w z,t where A is the amplitude of TME optimum oscillation and P is the period (the number of generations) of the TME cycle. The equation (21) allows the TME optimum to change periodically from 0 to A and then to 0 in P generations. In all simulations of this scenario P is fixed at 360 generations.

Stable (i.e. constant) optimum
We assume the microenvironment optimum is constant with z 0 opt :

Sudden change in the optimum
We assume the TME optimum suddenly changes from z 0 opt to z 1 opt : . CC-BY-NC-ND 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/377226 doi: bioRxiv preprint first posted online Jul. 25, 2018; z opt t ( ) = z 0 opt = z 0 opt ,...,0 ( ) , (25) z opt t ( ) = z 1 opt = z 1 opt ,...,0 In this scenario, the cancer cells evolve under constant stabilizing selection in the TME with optimum z 0 opt , and after a sudden change the cancer cells then evolve under another constant stabilizing selection with a new optimum z 1 opt (the optimum of the first trait changes from z 0 opt to z 1 opt ).

Simulation of cancer evolution in the three-dimensional space
Our intention was to identify key cancer evolutionary patterns in 3D and investigate how these are influenced by model parameters, such as the number of phenotypic traits (the complexity of cancer cells), changing microenvironment (properties of the TME), selection correlation, population size, mutation rate, initial phenotype/fitness of cancer cells (distance to optimum), initial population size and chromosome instability affect 3D cancer adaptation (tumorigenesis) in an ecological and evolutionary framework.
We simulate the cancer cell population in 3D space in one of the six environments described above (Supplementary Figure S1). The simulation is fully individual-based, following simple growth mechanics 10,60 . In brief, the offspring of newly divided cells stochastically look for available positions in a 3D lattice space. The population evolves in a discrete and non-overlapping manner. The empty 3D tumour space represents the tumour microenvironment that may change, following the dynamics described above, and the position in the 3D space does not affect the TME or selection.
Simulations are performed under particular TME change patterns due to either natural cause or anti-cancer therapies as described in equations (12)- (28). The number of traits is fixed in each simulation at n = 1 to n = 8 . Here we assume cancer cells require at least n . CC-BY-NC-ND 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/377226 doi: bioRxiv preprint first posted online Jul. 25, 2018; one adaptive trait for survival (e.g. anti-apoptotic or metabolism related traits) as also inspired by microbial adaptive evolution 61 and mammals with very low cancer incidencies 62,63 . The initial population has K neoplastic cells, and grows from the centre of the specified 3D coordinates in the 3D square lattice (it has 26-cell cubic neighbours, also called a 3D Moore neighbourhood) until it reaches the defined maximum population size or 3D tumour space. The initial K cells have an adjustable phenotype z 0 .
This determines the initial fitness in the fitness landscape, which has the optimum at the origin and an adjustable shape defined by S in equation (9) (illustrated in Supplementary Figure S2). The initial phenotype can also be calculated by equation (9) when the predefined fitness value is known. For generality we consider that the adjustable phenotype of initial K cells is determined by either classical driver mutations or by phenotypic plasticity.
The cancer cells have c sets of chromosomes (ploidy) and reproduce asexually with chromosome instability (CIN) (Supplementary Figure S1). To model CIN for simplicity we define a rate, r c , to vary the copy of paternal and maternal chromosome passed down to daughter cells. Each individual cancer cell is represented by L diploid loci (for computational efficiency without loss of generality c = 2 , L = 5 , unless stated otherwise ) with additive effect on the n-dimensional phenotype z .
A multifurcating tree tracks the genealogies of the alleles generated at these loci. When only one branch of an allele tree survives an adaptive step is recorded as per standard terminology 48 . The per locus mutation rate, u = 4 ×10 −5 , is used for per genome replication 64 . Due to high mutation rate and large population size multiple mutations may be generated in each generation. So there is clonal interference in our simulations. (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/377226 doi: bioRxiv preprint first posted online Jul. 25, 2018; to set the doubling time as 24 hours is to allow simulations of more generation time to gain more insights into the underlying evolutionary process (e.g., one day doubling time allows roughly 37000 generations for 100 years, however, 5     to their fitness cut-off value, which results in two colours for the 3D snapshots: red (large than or equal to cut-off value indicating higher fitness) and green (small than cutoff value indicating lower fitness), respectively. The cut-off value for fitness is 0.75 for (a, d). The cut-0ff value for fitness is 0.9 for (b, e). The cut-off value for fitness is 0.99 for (c,  a, the mean population fitness of the evolving cancer cells is sampled at various TME change rates . The dashed line represents mean fitness 0.5. When mean population fitness reaches this value it is more likely to be extinct. The line represents a simple linear regression fit. b, the mean selection coefficients are sampled at various TME . CC-BY-NC-ND 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/377226 doi: bioRxiv preprint first posted online Jul. 25, 2018; change rates. The line represents a smooth curve. c, the mean selection coefficients are plotted against the TME change rates. d, the mean population fitness is plotted against the number of traits of cancer cells. e, the mean selection coefficients are plotted against the number of traits of cancer cells. The TME rates used are v 1 = 0.5 , v 1 = 0.05 ,  Simulations are performed with different initial conditions regarding population size, fitness and TME change rates. Three conditions for initial population size are used, N = 10 (coloured blue), N = 10 4 (coloured red) and N = 10 7 (coloured green), respectively. Three conditions for initial fitness are also used, w = 0.1 (a-d), w = 0.5 (e-h) and w = 0.9 (i-l). Populations also evolve under four different TME change rates:  Example 3D snapshots are taken sequentially after sudden change of the TME optimum due to treatments (a-o). Four different treatments maintaining four different levels of TME optimum are shown: z 1 opt = 5 (a-e), z 1 opt = 6 (f-j), z 1 opt = 7 (k-o) and z 1 opt = 8 (no 3D snapshots are taken due to immediate population extinction after treatment), respectively. For understanding clonality of cancer cells after treatments, e, j and o show the clonal expansion of cancer cells after treatments for z 1 opt = 5 , z 1 opt = 6 and z 1 opt = 7 , respectively. In order to track the precise fitness status of the population the sample is taken for every generation. The population fitness plotted against generation time is . CC-BY-NC-ND 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/377226 doi: bioRxiv preprint first posted online Jul. 25, 2018; summarized in p for each treatment, and the dashed line indicates when the treatment starts (after generation 100).
. CC-BY-NC-ND 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/377226 doi: bioRxiv preprint first posted online Jul. 25, 2018; . CC-BY-NC-ND 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/377226 doi: bioRxiv preprint first posted online Jul. 25, 2018; . CC-BY-NC-ND 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.