## Abstract

Cancer dynamics are an evolutionary game between cellular phenotypes. A typical assumption in this modelling paradigm is that the probability of a given phenotypic strategy interacting with another depends exclusively on the abundance of those strategies without regard for local neighbourhood structure. We address this limitation by using the Ohtsuki–Nowak transform to introduce spatial structure to the go versus grow game. We show that spatial structure can promote the invasive (go) strategy. By considering the change in neighbourhood size at a static boundary—such as a blood vessel, organ capsule or basement membrane—we show an edge effect that allows a tumour without invasive phenotypes in the bulk to have a polyclonal boundary with invasive cells. We present an example of this promotion of invasive (epithelial–mesenchymal transition-positive) cells in a metastatic colony of prostate adenocarcinoma in bone marrow. Our results caution that pathologic analyses that do not distinguish between cells in the bulk and cells at a static edge of a tumour can underestimate the number of invasive cells. Although we concentrate on applications in mathematical oncology, we expect our approach to extend to other evolutionary game models where interaction neighbourhoods change at fixed system boundaries.

## 1. Introduction

The importance of heterogeneity within tumours is gaining ground as one of the most important factors in the laboratory and clinic alike [1]. This heterogeneity exists at multiple scales, each of which presents its own unique set of challenges. One form of phenotypic heterogeneity that has been widely studied was first described by Giese *et al.* [2], when they showed that in gliomas migration and proliferation were mutually exclusive processes; motile cells are unable to proliferate while they are moving and proliferating cells are unable to move while they divide. This has been termed the *go or grow* dichotomy [3]. A proliferative or autonomous growth (grow) cell might switch to a motile or invasive (go) cell either through mutation, metabolic stress [4], undergoing the epithelial–mesenchymal transition (EMT) [5], or some other mechanism. EMT is generally characterized by the loss of cell–cell adhesion and a gain in motility and invasiveness in tumour cells and is often the first step in the metastatic cascade. The opposite, or mesenchymal–epithelial transition, occurs when the cells arrive at a distant site and shift their invasive phenotype to one of more aggressive clonal growth [6]. This combination of transitions is one of the hallmarks of several carcinomas such as prostate [7,8], breast [9] and other ductal cancers, where pre-invasive neoplasms are constrained by architectural boundaries (edges) such as the duct wall or basement membrane. Our goal in this work is to highlight the important and overlooked role of such edges in the evolutionary dynamics of the competition between go and grow cells.

Evolutionary game theory (EGT) is a mathematical approach to modelling frequency-dependent selection, where players interact strategically not by choosing from a set of strategies but instead by using a fixed strategy determined by their phenotype. Given the evolutionary nature of cancer [10,11], EGT has been applied to study how the interactions between different types of cells in a polyclonal tumour could drive the dynamics of a given cancer [12]. In its first application to oncology [13,14], EGT was used to analyse the circumstances that lead to coexistence of two phenotypes. Subsequent research [15] extended this idea to interactions between three players in the angiogenesis problem. Gatenby and Vincent adopted a game theory model heavily influenced by population dynamics to investigate the influence of the tumour–host interface in colorectal carcinogenesis [16,17] and suggested therapeutic strategies [18]. Our own work, as well as that of others, has shown that EGT can be used to study the conditions that select for more aggressive tumour phenotypes in gliomas [19,20], colorectal cancer [16,17], multiple myeloma [21] and prostate cancer [22]. Furthermore, EGT has been used to investigate the impact of treatment on cancer progression [23–25]. For an in depth overview of the game-theoretical approach to cancer, see Basanta & Deutsch [12].

In this study, we introduce spatial structure into the canonical ‘go versus grow’ game [26,27] in which proliferation and motility compete within a tumour. We use our direct spatial approximation to consider a familiar scenario for conservation biology: the edge effect of an ecological system (for example, a forest in landscape ecology) at a static boundary [28–30]. In tumour progression, this is analogous to a cancer cell surrounded entirely by other cancer cells as opposed to constrained by a physical boundary, such as a basement membrane, organ capsule, fibrous capsule, or blood vessel (figure 1). Note that we do not consider the unconstrained, dynamically growing edge of the tumour.

The static boundaries studied in this article are exciting in and of themselves, as the evolutionary dynamics that occur at static boundaries govern progression past key cancer stages: the change from *in situ* to invasive; locally contained to regional advanced growth; and the dramatic shift from local to metastatic disease. The former situation occurs early in the progression of most epithelial tumours, and it is commonly believed that it is at this point when the Warburg effect occurs, pushing cells towards the glycolytic (acid producing) phenotype which promotes invasion and motility—the so-called ‘acid-mediated invasion’ hypothesis [19,31–34]. The latter situation occurs when tumours are up against blood vessels and is likely the first opportunity for haematogenous dissemination, the first step in the metastatic process [35]. Here, we show a striking change in the evolutionary game dynamics from the tumour bulk to the tumour's physical boundary (figure 1). This study represents, to our knowledge, the first attempt to understand the effects of changing neighbourhood structure on evolutionary game dynamics of tumours.

## 2. Material and methods

### 2.1. Inviscid game for motility

To mathematically model the *go or grow* dichotomy [3], we considered the situation in which the tumour is made of a population of rapidly proliferating cells capable of autonomous growth (AG), along with a subpopulation which arises by a mutation or phenotypic change which confers motility/invasiveness (INV) to tumour cells. The game has two parameters: *c* represents the direct and indirect costs of motility incurred by cells with the INV phenotype resulting from the reduced proliferation rate of motile/invasive cells [3] and *b* is the maximum fitness a tumour cell will have under ideal circumstances when it does not have to share space or nutrients with other cells. As the units of measure are arbitrary, the ratio *c*/*b* can be considered alone to determine game-theoretic dynamics of the proportion of cell types. With INV as strategy 1 and AG as strategy 2, the game's payoff matrix is
2.1

To understand the inviscid model [26], imagine two cells meeting at random in a resource spot; for an inviscid population this probability depends only on the cells' relative abundance but for the structured populations considered later, the probability of meeting again will be higher. If both cells are INV (motile), then one cell stays in the resource spot (i.e. a location that contains oxygen, glucose or other growth factors) and consumes the resources *b*, and the other pays a cost, *c*, to move and find a new empty site where it can then find resources, for a payoff *b*. If after a move a motile cell only finds a new empty spot with probability *r*, then the expected move payoff is *rb* − *c*. This can be captured as an indirect cost by adjusting the cost to be *c*′ = *c* + (1 − *r*)*b* without introducing an extra parameter. As the cell that has to move is chosen randomly from the two, the expected payoff for each cell is the average of the no-move (*b*) and move (*b* − *c*) payoffs. On the other hand, if an INV cell meets an AG cell then the INV cell will move, incurring the cost, *c*, before receiving the benefit, *b*, for a total fitness (*b* − *c*). The AG cell, however, will stay and consume all the resources (*b*). Finally, if two non-motile cells (AG) are in the same resource spot, then they simply share the resources, for a payoff of *b*/2.

Note that we are not explicitly modelling the mechanics of motility, but just their effects on cell fitness. Thus, our approach applies equally well to any other interaction with the same payoff structure. A more explicit account of the biomechanics of motility could inform the specific values of *b* and *c*, but the exact values of these parameters are not essential for our analysis and are beyond the scope of this study. We also expect the specific values of the payoffs to vary between tissues within a patient and across patients.

### 2.2. Direct approximation of spatial structure

A standard assumption in EGT is a perfectly mixed (inviscid) population, in which every cell in the population is equally likely to interact with every other cell [36]. This may be a reasonable assumption in some cases, but in most solid tumours (or any other situation being modelled) in which spatial structure is important, the validity of this assumption should be questioned [37]. The current solution to this is to map analytic EGT cancer models onto a lattice and run *in silico* experiments to simulate the resulting cellular automaton [26,38,39]. In such cases, the choices of the specific microdynamics to simulate are arbitrary and often left up to convention and the modeller's imagination, since direct empirical mechanisms at such a precise level are often unknown. Further, explicitly solving complex spatial structures is currently outside of existing mathematical tools, so computational approaches have to sacrifice the analytical power and full theoretical understanding of pure EGT approaches.

On occasion, computational modellers restore some analytical power by making mathematical approximations of the simulation that are already approximations of, or guesses at, the tumour's real spatial structure. To avoid this double approximation and to analytically model how spatial structure affects evolutionary games in the limit of large populations and weak selection, Ohtsuki & Nowak [40] derived a simple rule for taking a more direct first-order approximation of any spatial structure. This approach is based on the technique of pair-approximation [41–43] and is exact only for Bethe lattices (infinite trees of constant degree), but is highly accurate for any static structure where higher order terms, like the correlations between neighbours of neighbours, are negligible. For example, Ohtsuki and Nowak concentrated on the application of their tool to *k*-regular random graphs, which are locally tree-like and have negligible second-order and higher terms, but the transform can be used more broadly. While real spatially structured biological populations, such as solid tumours, can have non-negligible higher order interactions, Ohtsuki and Nowak's first-order approximation is an improvement over the common inviscid assumption that still allows us to explore a completely analytic model.

Given a game matrix *A*, one can compute the Ohtsuki–Nowak (ON) transform *A*′ = ON_{k}(*A*) and then recover the dynamics of the spatially structured game *A* by simply looking at the inviscid replicator equation of *A*′. Here, we present the transform in a form that stresses its important qualitative aspects:
2.2where ** 1** is the all ones vector and

**is the diagonal of the game matrix**

*Δ**A*= [

*a*], i.e.

_{ij}

*Δ**=*

_{i}*a*; thus,

_{ii}

*Δ*

*1*^{T}(

*1*

*Δ*^{T}) is a matrix with diagonal elements repeated to fill each row (column). The first summand is the original payoff matrix. The second summand accounts for the more common same-strain interactions that are a consequence of local dispersal. The type of perturbation in the third summand was shown by Hilbe [44] to result from finite sampling of interaction partners. The summands are not arbitrary, and emerge as a whole from the pair-approximation technique. Our rationale for grouping the summands in this particular form is to help build intuition for equation (2.2). Effects of indirect interaction from more distant cells could be introduced as further corrections to equation (2.2), but would require more assumptions about the tumour microstructure [43]. Because local interactions between neighbouring cells are empirically well studied [45], the focus of the current study is incorporating only this first-order structure into an EGT model of cancer dynamics.

## 3. Results

Whenever *b* > *c*, the game in equation (2.1) is a social dilemma (like Prisoner's Dilemma or Hawk–Dove game; for a classification, see [46]) with invasive cells as the cooperators, and AG cells as defectors. Rules of thumb from EGT [47] suggest that cooperators benefit from the structure of small interaction neighbourhoods, in agreement with our biological intuition that, in this game, having the ability for conditional motion is of more use in a more constrained and viscous environment than in one where all cells are already stochastically moving around and interacting at random. We look at this formally by explicitly considering spatial effects on the previously inviscid model. Applying the transform from equation (2.2) to the game in equation (2.1) yields
3.1

This game has three qualitatively different regimes that depend on the value of *c*/*b* and *k*:

(1) If (

*k*+ 1)/(*k*^{2}+ 1) ≥*c*/*b*, then there is a single stable fixed point with all cells invasive. All polyclonal tumours evolve towards this fixed point. For inviscid populations (*k*→ ∞), this condition is satisfied only if motility is cost-free (*c*= 0) and hence the possibility of an all invasive tumour was not noted in previous non-spatial analysis [26].(2) If (

*k*+ 1)/(*k*^{2}+ 1) <*c*/*b*< (*k*+ 1)/(2*k*+ 1), then the game has Hawk–Dove dynamics, and there is a stable polyclonal equilibrium with a proportion*p*of INV cells: 3.2 All polyclonal populations will converge towards this proportion of INV cells. In the unstructured limit as*k*→ ∞, we have perfect agreement with our previous results [26] and recover the condition*c*/*b*≤ 1/2 that was assumed for the inviscid equilibrium to exist and the exact numeric value of (*b*− 2*c*)/(*b*−*c*) for the equilibrium proportion of INV agents. For any finite*k*, however, the proportion of invasive cells is strictly higher.(3) If (

*k*+ 1)/(2*k*+ 1) ≤*c*/*b*, then the game has Prisoner's Dilemma dynamics and any polyclonal population converges towards all AG and the tumour remains non-invasive.

These three regimes are plotted in figure 2. When 1/(*k* − 2) = 0, we have the inviscid game and for 1/(*k* − 2) = 1 we have the most structured regime possible with small neighbourhoods (*k* = 3). The red region corresponds to a completely INV tumour, the yellow to a polyclonal tumour and the green to all AG. As we make the tumour more structured and reduce the number of neighbours, it becomes easier for the INV cells to be expressed in the tumour.

For example, consider the case where *c*/*b* = 0.53, and as in figure 1, in the tumour away from boundary (the turquoise cell) there are *k* = 8 neighbours and the dynamics favour all AG, so no INV cells will be present at equilibrium. For an example of replicator dynamics of this condition, see the green line of the right inset. If the solid tumour is pressed up against a static boundary then cells at the edge have fewer neighbours (e.g. *k* = 5, the pink cell) and the dynamics at the boundary favour a polyclonal population with about 8% of the cells INV. For an example of replicator dynamics for this condition, see the yellow line of the right inset. Note that the tumours represented by the green (*k* = 8) and yellow (*k* = 5) lines of the right inset have the same *c*/*b* = 0.53 and initial proportion of invasive cells *p*_{0} = 0.04 and yet the invasive phenotype is pushed to extinction in the tumour bulk (*k* = 8), but stabilizes near potentially dangerous level of invasive cells (*p* = 0.08) at the tumour edge (*k* = 5). Similar higher selection for invasiveness at the static edge is present for the more competitive environment of *c*/*b* = 0.23 in the left inset, but in this case we started the example replicator dynamics at *p*_{0} = 0.93. In this case, we have a polyclonal tumour bulk (*k* = 8) with *p* = 0.86 and convergence towards all invasive phenotypes at the static edge (*k* = 5). Of course, the specific parameter values above are for illustrative purposes. Actual change in *k* (just like *c*/*b*) will be experiment and geometry dependent—for example, we expect a more drastic decrease in *k* for a convex rather than a concave boundary.

The main message is that the edge effect can cause a polyclonal boundary in a tumour with a homogeneous all AG body. And while not yet shown to be universally applicable, this is in qualitative agreement with the experience of pathologists, such as the typical staining for motility (EMT) in figure 3, where we see the bone marrow space entirely invaded by carcinoma cells (pancytokeratin; figure 3*a*), but a stark difference in motile phenotype as illustrated by the nuclear staining by SLUG, which is only present in high concentration in cells against the static boundary of the bony trabecula. As with all interesting questions in biology, there can be many other causes involved in such specific experimental or clinical systems, and our result is focused on the idealized change in local architecture common to many of them. The effect to which we are drawing attention is different from and can occur in addition to other system-specific mechanisms or heterogeneity due to more traditional sources like gradients in metabolites and drugs.

## 4. Discussion

A standard assumption in EGT is that all players interact with all other players: the population is inviscid. There are a number of biological scenarios in which such an assumption could be misleading, and we consider such a scenario in the form of several key aspects of solid tumour progression. The role of spatial heterogeneity has been explored before in mathematical models studying evolutionary processes in cancer [48–50]. Those models show that different environments produce different selective pressures and that phenotypic heterogeneity results from, and drives, the spatial one. Our approach tackles a different question related to the nature of physical edges on cancer evolution for which we have applied the Ohtsuki–Nowak transform [40] to the standard go–grow game of mathematical oncology [26]. We have shown a quantitative effect: spatially structured tumours promote the invasive phenotype compared to inviscid tumours with the same ratio of cost of motility to benefit of resources *c*/*b*.

We also considered the decrease in neighbourhood size experienced by cells at the static boundary of a tumour, compared to cells within the tumour bulk. This could represent a number of feasible and very relevant scenarios, including a tumour against a blood vessel, organ capsule or fibrous capsule; or an *in situ* neoplasm at the basement membrane. We have shown that this change in neighbourhood for tumour cells can, independently from any other parameters, significantly affect the evolutionary stable strategies: in this case, a promotion of the INV phenotype. The edge effect allows a tumour that internally has no invasive phenotypes expressed to have a polyclonal boundary with both invasive and non-invasive cells, a scenario which is known to appear in several cancers, most notably prostate [51] in the form of peri-neural invasion, and the recently described perivascular invasion in melanoma [52] (previously called angiotropism). In each of these two scenarios, tumour cells express the motile phenotype, but only when against the physical structure in question. Another recent result suggesting the importance of effects of changes in local architecture is the work of Belmonte-Beitia *et al.* [53], in which the authors show that the dynamics of motility change drastically at the grey–white interface in glioblastoma in ways that cannot be predicted by simply comparing the behaviour in each zone individually.

The results of our mathematical model could have significant translational implications. Genetic heterogeneity has recently become recognized as the rule in cancer [54], but for as long as physicians have had microscopes, we have realized that spatial organizational heterogeneity was an equally defining factor. The Gleason score [55] is a classic example of greater heterogeneity predicting worse survival in prostate cancer. We have shown that a specific change in neighbourhood size at a static boundary can dramatically alter game dynamics, and select for novel phenotypes, and gives a rationale for working to understand within-tumour differences, not just at relatively long length scales [56], but also at architecturally different locations at short length scales, which can be done with the advent of single-cell technologies and laser capture microdissection [57].

We have shown that the local spatial structure of a tumour can strongly affect the evolutionary pressures on its constituent cells, even if all other factors are held constant. This can add yet another source of sampling bias to tissue biopsies and suggests that the architectural, not just the molecular, context is important. For instance, consider an idealized fine-needle aspiration biopsy [58], assuming the standard 0.7 mm needle samples a perfect column around 20 cells in diameter of tumour cells right next to a critical boundary such as a capillary. In our running example of *c*/*b* = 0.53 (figures 1 and 2), this would result in the sample containing only about 0.4% invasive cells since 19 out of every 20 cells are not at the boundary. This is below the detection levels of the state-of-the-art medical practice [59,60]. However, among the critical 1 out of every 20 cells at the boundary, a dangerous proportion of 8% would display the invasive phenotype. Thus, an oncologist performing a diagnostic fine-needle aspiration biopsy could be led to think that a tumour poses a low risk for invasion/metastasis, because the technique destroys the local structure of the tumour and mixes the cells at the critical tumour boundary with the (in this case) irrelevant tumour bulk.

Our model was motivated by the study of cancer, but the spatial edge effects in games of interacting players that we investigate could represent any number of other scenarios. While we have focused on the specifics of metastasis and cancer invasion, this method could yield insights into many other interesting problems ranging from ecology to medicine, and highlights the importance of neighbourhood geometry when studying the evolutionary dynamics of competing biological agents.

## Authors' contributions

A.K., J.G.S. and D.B. designed the study, interpreted the results and wrote the manuscript. A.K. and J.G.S. carried out the mathematical modelling. All authors gave final approval for publication.

## Competing interests

We declare we have no competing interests.

## Funding

D.B. would like to acknowledge U01CA151924–01A1 for financial support. J.G.S. would like to thank the NIH for its generosity in providing a Loan Repayment Grant.

## Acknowledgements

We would like to thank Colm Morrissey from the University of Washington and Marilyn Bui from the Moffitt Cancer Center for illustrative images used in this paper, and David Liao for extensive helpful feedback.

- Received February 19, 2015.
- Accepted May 12, 2015.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.