## Abstract

Cancer evolution and progression are shaped by cellular interactions and Darwinian selection. Evolutionary game theory incorporates both of these principles, and has been proposed as a framework to understand tumour cell population dynamics. A cornerstone of evolutionary dynamics is the replicator equation, which describes changes in the relative abundance of different cell types, and is able to predict evolutionary equilibria. Typically, the replicator equation focuses on differences in relative fitness. We here show that this framework might not be sufficient under all circumstances, as it neglects important aspects of population growth. Standard replicator dynamics might miss critical differences in the time it takes to reach an equilibrium, as this time also depends on cellular turnover in growing but bounded populations. As the system reaches a stable manifold, the time to reach equilibrium depends on cellular death and birth rates. These rates shape the time scales, in particular, in coevolutionary dynamics of growth factor producers and free-riders. Replicator dynamics might be an appropriate framework only when birth and death rates are of similar magnitude. Otherwise, population growth effects cannot be neglected when predicting the time to reach an equilibrium, and cell-type-specific rates have to be accounted for explicitly.

## 1. Introduction

The theory of games was devised by von Neumann & Morgenstern [1], and according to Aumann [2], game theory is an ‘interactive decision theory’, where an agent's best strategy depends on her expectations on the actions chosen by other agents, and vice versa. As a result, ‘the outcomes in question might have been intended by none of the agents’ [3]. To rank and order strategies, and to optimize individual payoffs, different systems to systematically identify equilibria have been defined. Most famously, the Nash equilibrium is a set of strategies such that no single agent can improve by switching to another strategy [4]. This concept includes mixed equilibria, which describe probability distributions over strategies. Such equilibrium concepts in game theory cover various kinds of patterns of play, i.e. simultaneous, non-simultaneous and asymmetric strategies [5]. This rich and complex framework allows for a wide application of game theory beyond economics, famously in ecology and evolution [6]. In biological context, and especially in evolutionary game theory, the focus has been on simultaneous and symmetric strategic interactions in evolving populations [7].

Evolutionary game theory replaces the idea of choice and rationality by concepts of reproduction and selection in a population of evolving individuals [8] and was conceived to study animal conflict [9]. Behavioural phenotypes are hardwired to heritable genotypes. Without the possibility of spontaneous mutation events, offspring carry the parent strategy. Evolutionary games have also been used extensively to study learning and pairwise comparison-based changes in strategy abundance in populations of potentially erroneous players [10–12].

Selection in evolutionary games is based on the assumption that payoff translates into Darwinian fitness, which is a measure for an individual's contribution to the pool of offspring in the future. Complex deterministic dynamical systems arise when one considers very large populations of reproducing individuals. The most prominent example for such a system is the replicator equation [13], which focuses on the relative abundance of each strategy. The replicator equation does not model population growth specifically, but rather describes changes in relative abundances. Existence and stability of fixed points in these dynamical systems depend on the payoffs [14], and on the choice of fitness function [15]. In the study of animal behaviour, the precise measurements of payoffs, as observed from individuals' behaviours, is difficult. Milinski *et al.* [16] determined all but one payoff parameter precisely, in order to observe tit-for-tat strategies in repeated Prisoner's Dilemma games in fish. Kerr *et al.* showed that *Escherichia coli* bacteria can be observed to evolve according to rock–paper–scissors type of interactions, if cellular dispersal is minimal. A recent expansion of interesting theoretical considerations that apply evolutionary games to biology [17] occurred because of the ability to assess many problems in ecological and evolutionary population dynamics at least in qualitative terms, i.e. by predicting and ranking evolutionary equilibria, how population-wide coexistence can emerge from apparent individual conflict, or how fast transitions between equilibria occur.

Tumour cell populations, including cells of the tumour microenvironment, are part of a complex ecosystem [18], which can have consequences for therapeutic outcomes [19]. At the same time, it has been more widely recognized that Darwinian selection plays a key role in cancer [20]. Given the appreciated amount of both genetic and phenotypic heterogeneity in tumour cell populations [21], evolutionary games have become more widely used as a means to theoretically model tumour evolution, especially after tumour initiation [22]. Prominent examples of recent applications of replicator equations in cancer are concerned with the avoidance of the tragedy of the commons, where a sub-population of tumour cells produces a ‘tumour public good’ in form of an insulin-like growth factor [23], in form of glycolytic acid and vascular endothelial growth factor [24], or modelling the dynamic equilibrium between lactate respiration and glycolysis in tumour cells [25]. Such non-autonomous effects between tumour cells had been proposed some time ago [26], and non-cell-autonomous growth rates were recently measured empirically [27]. Similar findings and future challenges in this field have been summarized by Tabassum & Polyak [18].

We here focus on the time it takes to reach an equilibrium in different approaches to model deterministic evolutionary game dynamics. In particular, we focus on differences between logistic growth and the replicator dynamics. We show that the time to get arbitrarily close to an equilibrium, which we here call the *ɛ*-fixation time, might critically depend on the underlying cellular birth and death rates. We focus on two coevolving tumour cell populations, and present a discussion of the dynamics between growth factor producers *C*_{1} and free-riding non-producers *C*_{2}. In the simplest setting, we can assume that these closely related cell types experience population doubling rate *α*, and that the tumour public good, produced by *C*_{1} cells, has a linear positive effect on cellular birth rates in form the additive benefit that scales with the doubling rate *βα*, but bears a production cost *κ*. The respective game can be recast in the payoff matrix
1.1We assume that the linear benefit of the public good arises through growth factor diffusion that occurs on a time scale much faster than the average times between cell divisions. In a well-mixed population with fraction *u* of *C*_{1} cells, the fitness functions of this simple game are then be given by
1.2Our analysis in this paper is based on cell-type-specific doubling rates, and in the case of logistic growth, also on the apoptotic rates. We are interested in the question of when replicator dynamics, that typically only models changes in relative abundance as a result of fitness differences *f*_{1}(*u*) − *f*_{2}(*u*), predicts similar *ɛ*-fixation times as a logistic growth dynamics, and when this is not the case. The main idea is that the replicator dynamics neglects apoptotic rates, but that these rates in turn influence the time to reach an equilibrium in a co-growing and coevolving heterogeneous cell population.

## 2. Methods

In this section, we introduce our model of bounded frequency-dependent growth. We define our basic deterministic framework of two co-growing cancer cell populations, derive dynamic equations for the fraction of one clone and the total size of the population, and then derive an expression for the stable manifold of the system.

### 2.1. Logistic growth model

The population is assumed to consist of two types, and we denote their absolute numbers by *x*_{1} and *x*_{2}. The carrying capacity is denoted by *K*, which we consider to be a constant. It is possible to model it as a function of the strategies present in the population [28,29]. The growth rate of each type is assumed to depend on the fraction of type 1 cells *u* = *x*_{1}/(*x*_{1} + *x*_{2}) according to growth functions *f*_{1}(*u*) for type 1 and *f*_{2}(*u*) for type 2. Lastly, cells of both types die at a constant rate *μ*. Taken together this implies that we get the following system of coupled logistic equations that describe co-growth and coevolution of the two cell types:
2.1defined for . In the following analysis, we first assume *μ*_{1,2} = *μ* and *f*_{1,2}(*u*) > *μ* for *u* ∈ [0, 1], i.e. the net growth rate of both cells types will always be positive. In the second part of the discussion, we will relax the assumption of equal rates and turn to the more general case of *α*_{1} ≠ *α*_{2}, *μ*_{1} ≠ *μ*_{2}, as we analyse the system implementing previously measured cellular rates of proliferation and apoptosis. Note that the logistic growth model emerges from a spatial setting that includes cell movement if cell migration occurs on a much faster time scale compared to cell division. It has been shown that in this case spatial correlations are negligible and the population dynamics can be described using a logistic growth equation [30]. In this parameter regime, it is also justified to assume that interactions that influence the rate of cell division become independent of specific local configurations, and depend solely on the frequency of different cell types.

### 2.2. Analysis

To simplify the analysis of the system (2.1), we apply the following change of variables:
2.2where *u* is the fraction of type 1 cells and *s* is the total population size. By differentiating *u* and *s* with respect to time we obtain the following system of ODEs:
2.3defined on *u* ∈ [0, 1] and . We note that in the case when *s* is small compared with the carrying capacity *K*, such that *s*/*K* ≈ 0 the system reduces to
2.4and we see that the equation for *u* is independent of the population size *s* and that *u* changes according to the standard replicator equation [13,14]. We will now proceed to a more general analysis of our model.

#### 2.2.1. Fixed points

By solving the equations
2.5we see that for all growth functions *f*_{1} and *f*_{2} the system has the following set of fixed points on the boundary (see appendix A for details):

(1) (

*u*_{1},*s*_{1}) = (0, 0) with corresponding eigenvalues*λ*_{1}=*f*_{1}(0) −*f*_{2}(0) and*λ*_{2}=*f*_{2}(0) −*μ*> 0, which is unconditionally unstable,(2) (

*u*_{2},*s*_{2}) = (1, 0) with corresponding eigenvalues*λ*_{1}=*f*_{1}(1) −*f*_{2}(1) and*λ*_{2}=*f*_{1}(1) −*μ*> 0, which is unconditionally unstable,(3) (

*u*_{3},*s*_{3}) = (0,*K*(1 −*μ*/*f*_{2}(0)) with corresponding eigenvalues*λ*_{1}= (*μ*/*f*_{2}(0))(*f*_{1}(0) −*f*_{2}(0)) and*λ*_{2}=*μ*−*f*_{2}(0) < 0, which is stable iff*f*_{1}(0) <*f*_{2}(0), and(4) (

*u*_{4},*s*_{4}) = (1,*K*(1 −*μ*/*f*_{1}(1)) with corresponding eigenvalues*λ*_{1}= (*μ*/*f*_{1}(1))(*f*_{2}(1) −*f*_{1}(1)) and*λ*_{2}=*μ*−*f*_{1}(1) < 0, which is stable iff*f*_{2}(1) <*f*_{1}(1).

Here, fixed points 1 and 2 are trivial in the sense that they correspond to a system void of cells. Fixed points 3 and 4 correspond to monoclonal populations and are stable if the resident type has a larger growth rate compared with the invading type.

If there are points *u** ∈ (0, 1) such that *f*_{1}(*u**) = *f*_{2}(*u**), then these give rise to fixed points (*u**, *K*(1 − *μ*/(*f*_{1}(*u**)*u** + *f*_{2}(*u**)(1 − *u**)))), which are stable if *f*′_{1}(*u**) − *f*′_{2}(*u**) < 0 (see appendix A for proof).

We note that the stability criteria for the non-trivial fixed points at *u* = 0 and 1, including potential internal fixed points, are identical to those of the two-type replicator equation with payoff functions *f*_{1} and *f*_{2}.

#### 2.2.2. Invariant manifold

We now focus our attention to the dynamics when the system is close to saturation (*s* ≈ *K*) with the aim of obtaining a simpler description of how the frequency *u*(*t*) changes in time. This can be achieved since the phase space contains a stable invariant manifold that connects all the non-trivial steady states. The invariant manifold is simply a curve *s* = *h*(*u*), which attracts the dynamics and once the system enters the manifold it will not leave it. This implies that the dynamics along the manifold is effectively one-dimensional, and can be captured with a single ODE for *u*(*t*).

If we write the invariant manifold as a function *s* = *h*(*u*), then, since it is invariant it must be tangent to the vector field (d*u*/d*t*, d*s*/d*t*) at every point. This implies the condition
2.6which is known as the manifold equation [14,31]. By substituting d*u*/d*t* and d*s*/d*t* from (2.3) and letting *s* = *h*(*u*), we obtain the following equation for *h*(*u*):

This equation is a nonlinear ordinary differential equation and in order to solve it we express *h*(*u*) as a series expansion in the death rate *μ*, which typically is a small parameter
2.7where *a*_{i}(*u*) are coefficients that depend on *u*. We insert this ansatz into equation (2.6) and equate powers of *μ* to solve for the *a*_{i}'s. We do this for *i* = 0, 1, 2, introduce , and get
Numerical comparison shows that the invariant manifold is closely approximated by the first two terms, and we therefore drop all higher order terms and approximate the invariant manifold with
2.8Note here that the complete solution would be more complicated, as can be seen from the fact that this first order expression does not solve the original manifold equation.

The dynamics along the invariant manifold are given by replacing *s* with in (2.3), and we get the following expression (to first order in *μ*):
2.9With the unusual prefactor that is inversely proportional to the total fitness of the population, , this equation for the frequency of type 1 cells is similar to the version of the replicator equation introduced my Maynard-Smith [32], and the one derived by Traulsen *et al.* [33] (if we disregard the demographic noise term). The difference compared to previous derivations is the factor *μ*, which implies that the rate of change of *u* along the invariant manifold is proportional to the death rate.

## 3. Results and discussion

It is often argued that prefactors to the replicator equation are irrelevant since the dynamic flow and fixed points remain unchanged. However, the time scale of selection leading to an equilibrium might be altered. In this section, we explore the difference between the standard replicator equation and the logistic model considered here. We examine this relationship in the context of a tumour public goods game, in which some cells produce a public good at a cost, rendering a benefit to *all* cells in the population.

### 3.1. Diffusing public goods game

Autocrine production of growth factors is a common feature of cancer cells, and has previously been modelled using evolutionary game theory [23,34]. Let us now consider two cell types that only differ in one aspect. Type 1 cells produce growth factor at a cost *κ*. Type 2 cells do not produce the growth factor and are termed free-riders. Otherwise, both cell types have the same growth rates, which are a linear function of growth factor availability. We assume that the growth factor production rate is given by *ρ* and that the growth factor is bound and internalized by both cell types at rate *δ*.

Two largely simplifying assumptions are that, first, we are describing a well-mixed system and that, second, the growth factor concentration *G* is assumed to be uniform in space. We rely on the first assumption for mathematical convenience, as otherwise we would have to resort to non-analytical, agent-based or hybrid modelling [35]. Secondly, additional growth factor provision was shown to be rapid and leading to high levels of tumour public good, provided that the respective genetic promoter was strong [27]. In a similar study by Cleary *et al.* [36], who studied Wnt1-based cooperative tumour evolutionary dynamics, aberrant expression of the cooperative signalling molecule was observed on a tumour wide scale. Thus, under these simplifying but productive assumptions, the growth factor dynamics obeys the equation
Further, we assume that the growth factor dynamics occur on a fast time scale compared to changes in *x*_{1} and *x*_{2}. This implies that
and we can solve for *G* to give
where *β* = *ρ*/*δ*. For simplicity, we first consider a linear effect of the growth factor on the rate of cell division, as well as equal proliferation and death rates, which results in the growth functions given by equations (1.2). In order for the growth rate to be larger than the death rate for all *u* we assume the inequality *α* − *κ* > *μ*. This choice of growth functions gives the following system of ODEs for the frequency of producers *u* and the total population size *s*:
3.1This system results from equations (2.3) and has two non-trivial steady states given by a monomorphic population of free-riders (0, 1 − *μ*/*α*), and a population consisting only of producers (1, 1 − *μ*/(*α*(1 + *β*) − *κ*)); see analysis following equations (2.5). The eigenvalues are
3.2and
3.3and hence the free-rider steady state is stable. For the other fixed point (producers dominate), we have
3.4and
3.5making it unstable. Figure 1*a* shows the phase space of the system, where the open circles indicate unstable steady states and the filled circle shows the location of the single stable steady state. We note that for almost all initial conditions the dynamics rapidly converge to the invariant manifold (2.8) which is approximately given by
3.6Once the system enters the invariant manifold the dynamics can be approximated by (2.9) which for the diffusing public goods game considered here are given by
3.7Thus, in order to assess the impact of cell death and turnover on selection, we compare our description of the public goods game (3.1) with the standard replicator equation
3.8Figure 1*b* shows a comparison between the solution of the logistic system (3.1) and the replicator equation (3.8) for the same initial condition *u*_{0} = 0.75 (*s*_{0} = 0.01*K*) and with a death rate of *μ* = 0.1 h^{−1}. Whereas the two solutions agree for small times (when *s*≪*K*), they start to diverge as soon as the solution to the logistic system enters the invariant manifold. The solution of the replicator equation quickly converges to the steady state *u* = 0, while the fraction of producers in the logistic case decreases approximately linearly with time.

To quantify the effect of the death rate *μ* on the rate of selection we measured the time it takes for the logistic system to approach a steady state. For a fixed initial condition (*u*_{0}, *s*_{0}) = (0.75, 0.01), we measured the time it took for the system to reach a small *ɛ* neighbourhood of the fixed point, i.e. |*u*(*t*) − *u**| ≤ *ɛ*, with *u** = 0 and *ɛ* = 0.01. We call this time the *ɛ*-fixation time. All other parameters were fixed at *α* = *β* = 1, *κ* = 0.1, *μ* = 0.1 and *K* = 1. The result is displayed in figure 1*c* and shows that the *ɛ*-fixation time scales as *μ*^{−1}. This implies that for small *μ* the time it takes the system to reach the steady state can be exceedingly long. It is worth noting that the *ɛ*-fixation time for the replicator equation can be obtained in the limit of , performed on the logistic system, implying a never-growing population, in which the death rate equals the average birth rate.

### 3.2. Time scales of *in vivo* and *in vitro* cellular expansions

Previous studies of ecological interactions in growing tumour cell populations have observed various forms of frequency-dependent effects. These effects have then been linked to the persistence of distinct cancer cell lines that provide growth enhancing public goods to the tumour, most notably in experimental work by Marusyk *et al.* [27]. There, it could be shown that a mixture of certain clones could not explain tumour outgrowth *in vivo* by simply using superposition of individual clonal birth and death rates. Rather, synergistic tumour-driving effects can emerge, pointing to more intricate, potentially frequency-dependent growth effects, based on direct or indirect clonal interactions [18]. For the purpose of illustration, we extracted individual clonal proliferation (*α*_{i}) and death rates (*μ*_{i}) from Marusyk *et al.* [27], in order to predict how these rates shape the dynamics. Out of 16 clonal cell lines, each distinctively expressing a different gene, we chose four clones to calculate baseline cellular birth and death rates. The four clones, derived from the breast cancer cell line MDA-MB-468, were LoxL3 (lysyl oxidase type 3 [37], linked to breast cancer invasion and metastasis), IL11 (interleukin 11, a member of the IL 6 family that plays a multifaceted role in leukemia and breast cancer [38]), and CCL5 (C-C motif ligand 5, a chemokine with emerging roles in immunotherapy [39]). The baseline cellular birth and death rates of these clones were calculated in the following way, based on *in vivo* growth experiments, originally performed in a mouse xenograft model (tumours formed by orthotopic trans-plantation into the mammary fat pads of immunodeficient Foxn1^{nu} (nu) mice [27]). For all four clones, it was established that tumours grew exponentially; from longitudinal measurements and associated cellularity calculations, the net cellular doubling rates were calculated (see Ext. Data fig. 3 and SI in [27], where exponential growth rates are given, which we transformed into doubling rates). For the four above mentioned clones, proliferation assays were also performed (Ext. Data fig. 1 in [27]). These BrdU staining experiments measure the fraction of cells in S-phase of the cell cycle, *χ*. S-phase duration *T*_{S} is highly conserved in mammary cells [40], known to be approximately 8 h long, *χ* serves as a direct estimate for the per cent of S-phase in relation to the whole cell cycle *T*, and thus the doubling rate, which we set to *α* = 1/*T*. Using the relation
3.9we calculated the monoclonal birth rates using
3.10Thus, given the net doubling rate *r* = *α* − *μ*, it is possible to estimate the death rate
3.11with *T*_{S} fixed to 8 h. Data for *r* and *χ* are given in appendix B. As for both *r* and *χ*, several independent measurements were performed, and we calculated distributions of *α* and *μ* for the three cell lines described above. We contrasted these distributions to *in vitro* distributions of cellular birth and death rates, adapted from [41] (fig. 3 therein), which are, notably, very similar to other *in vitro* values, e.g. reported for the PC-9 non-small cell lung cancer cell line [42] (figure 2*a*). In the *in vivo* tumour growth experiments, exponential growth was observed within the time frame of 50–80 days, at growth rates up to two population doublings per day (net growth rate) [27]. However, in most tumours the net growth rate was more moderate, and the actual cellular birth and death rates were at least of similar order in magnitude (*α*/*μ* ≈ 1). This stands in contrast with the birth–death rate ratios observed in cell cultures, where birth rates often exceed death rates by an order of magnitude (*α*/*μ* ≈ 10) [41–44].

As a notable difference to the previous section, here we assume both *α*_{1} ≠ *α*_{2} and *μ*_{1} ≠ *μ*_{2}. Thus, instead of (2.3), we now deal with the more general payoff structure
3.12and obtain the following ODEs for frequency of producer cells and total size of the system:
3.13and we seek to estimate the time it takes to reach a small *ɛ* neighbourhood of the equilibrium |*u*(*t*) − *u**| ≤ *ɛ*, shown in figure 2*b*. The combination of IL11 with one other cell line was chosen because it has been established that IL11 is a growth factor producer clone, which, at least in a first approximation, renders a linear fitness benefit [27]. We here make the additional assumption that IL11 cells carry a cost associated with growth factor production, and explore the extinction process of IL11 cells as they compete with either CCL5 or LoxL3 cells (figure 2).

We can calculate an estimate of the ‘time to fixation’ in the following way. Suppose the fraction of growth factor producers, *u*, is at a stable equilibrium, and that there are only two possible stable equilibria, *u** = 0 and *u** = 1. Then, the stationary solutions for the population size, *s**(*u**), will be
3.14and
3.15We now assume that the total population size remains at the stationary value, although it in fact changes (slightly) with *u*. This assumption can be thought of as a zeroth-order approximation in 1 − *s*/*K*, and it implies that near the stable manifold, the frequency *u* obeys the ODE
3.16which we can solve by inserting the approximations (3.14) and (3.15) into the ODE (3.16) and get the two solutions (for two different possible endpoints)
3.17and
3.18We now seek solutions of |*v*_{0,1}(*τ*) − *u**_{0,1}| ≤ *ɛ* for *τ* (with the equilibrium points *u**_{0} = 0, *u**_{1} = 1), and find the following relations that approximate the *ɛ*-fixation times:
3.19where *u*_{0} is the initial frequency. Note that here, we deviate from the notion of (average) fixation times in the strict stochastic sense [45], and replace the term by a threshold-based analytical approximation. Especially in a population that has reached the stable manifold, even a small fraction of remaining producer cells could still mean that there are as many cells as needed to warrant a mean field rather than a fully stochastic description.

For the *u* → 0, *s* → *K*(1 − *μ*_{2}/*α*_{2}) case, we can now compare our analytical approximations with the *ɛ*-fixation times of the full numerical solution in figure 2*c*, as a function of *κ*. Depending on the differences in clonal birth and death rates, the approximation exhibits qualitative differences. Equation (3.19) consistently overestimates the *ɛ*-fixation time if the death rate of the producer cells is lower than that of non-producers (IL11 with CCL5, *α*_{1} − *μ*_{1} < *α*_{2} − *μ*_{2}), but it underestimates the *ɛ*-fixation time if the net growth rate of the producer cells is higher than that of non-producers as long as the cost of growth factor production does not exceed a certain threshold (IL11 with LoxL3, *α*_{1} − *μ*_{1} ≈ *α*_{2} − *μ*_{2}). Hence, not only the cost of growth factor production factor influences the time to extinction of producer cells, but also the monoclonal net growth rate influences both the time to extinction of producers and the impact of an assumed cost associated with growth factor production. The approximations (3.19) are of ‘zero-order’ in changes in *s*. Yet, they are still able to reflect the basic fact that *ɛ*-fixation time can be heavily influenced by the cellular death rate of the abundant cell type. According to our rough approximation, the extinction time of producer cells (3.19) is both proportional to the ratio of birth to death rate of the non-producers, as well as inversely proportional to the birth rate difference. Surprisingly, in this approximation *τ*_{u →0} does not depend on the absorption or production rate of the growth factor, captured by *β*. Large differences in baseline birth rates extend growth factor producer extinction times. For larger values of *α*_{2}/*μ*_{2}, the extinction time is less sensitive to changes in the cost of growth factor production.

The two cellular death rates *μ*_{1} and *μ*_{2} have different effects on *ɛ*-fixation times. We used numerical solutions of the full system (3.13), in comparison to the replicator equation (3.8), to analyse variability of *ɛ*-fixation times (extinction of growth factor producer cells) under variable individual death rates. Thereby, we recover that higher total death rate speeds up the *ɛ*-fixation time across different initial conditions (figure 3*a*), and that the death rate of the ‘winner-clone’ plays a more important role (figure 3*b*): *μ*_{2} has a more pronounced impact on the *ɛ*-fixation time of non-producers. This might be connected to the fact that apoptosis-driven cell turnover of the nearly dominant cell type (i.e. the non-producer cells) governs the *ɛ*-fixation time. In accordance with this observation, the stable manifold is itself governed by the apoptotic rate of the dominant clone; compare to equation (2.8).

## 4. Summary and conclusion

We here have presented calculations that were concerned with the stability and time to reach a neighbourhood of equilibrium points in evolutionary game dynamics between two types of tumour cells. We focused on the dynamics of a tumour public good (tumour growth factor), in which we assumed linear fitness functions of growth factor producers and non-producers. The fitness function linearly depends on the relative abundance of growth factor producers, and production comes at a cost. We did not assume that the evolving population was at carrying capacity, as reflected in the logistic growth model. Thus, in general, population expansion and cellular birth as well as death rates are of importance for the time the system takes to equilibrate. The standard replicator equation typically rules out explicit death effects, and thus may not accommodate the impact of these death rates on the time to reach a population equilibrium.

The use of replicator equations and birth–death processes assume constant population size [7] or a population which is growing uniformly, for example, at an exponential rate [13]. These assumptions have led to a plethora of fruitful results in evolutionary game theory [46], e.g. to the ability to understand fixation and extinction times in evolutionary 2 × 2-games [47–50], multiplayer games [51], structured populations [52] or bi-stable allelic competition [53,54]. Evolutionary games have also been used to establish rules for equilibrium selection even in complex group-coordination games [55,56], in chemical game theory [57], and to map complex tumour dynamics [23–25,34,58–61]. However, the assumption that the population is either at constant size may be limiting, as also recently discussed by Li *et al.* [62] in the context of co-growing and coevolving bacterial species. Instead, the near-equilibrium population size and the time to reach equilibria are influenced directly by birth and death rates in the population.

We show that, for small differences between the birth and death rates, the eco-evolutionary dynamics of the mixture of two clones may be approximated by standard replicator dynamics. Analysis of previously established growth factor-dependent tumour dynamics of *in vivo* tumour growth showed that this parameter regime might indeed be biologically relevant (figure 2), even when the tumour population has not reached its carrying capacity. However, prominent examples of *in vitro* cell line expansions demonstrate that large differences between cellular death and birth rates might impact the dynamics in a different way [42–44], and in this case the replicator equation is a poor approximation of the eco-evolutionary dynamics. We used a logistic growth model that includes cell death. This system describes both co-growth, as well as coevolution of two tumour cell types. The choice of logistic growth is by no means unique, but a simple, first-order form of non-uniform growth.

We report two major findings. First, a first-order approximation in death rates allows estimation of the stable manifold, and reveals linear dependence on the apoptotic rate of the more abundant cell type. Second, this knowledge can be used to inform a zero-order approximation (in constant system size) of the time to get arbitrarily close to equilibrium (*ɛ*-fixation time), which reveals that indeed the cellular turnover of the dominant cell type near equilibrium governs the *ɛ*-fixation time as the system slowly moves along the stable manifold. This framework allowed us to examine the degree of the resulting variability in *ɛ*-fixation times based on previously measured *in vivo* tumour cell proliferation and death rates in the context of competition between producers and non-producers of a growth factor public good.

Various aspects of cancer cell population structure, such as cellular differentiation, localization or spatial heterogeneity, point to dynamic nonlinear size changes over time, especially during treatment [63–66], and treatment can shift the evolutionary game [67]. Furthermore, selection mechanisms that go beyond relative fitness differences play a role in our understanding of other biological and clinically relevant systems, such as the hematopoietic system [68,69]. Hence, future modelling efforts that seek to apply evolutionary game theory to explain complex cancer growth patterns need to precisely disentangle complex interaction patterns between cells from the overall growth kinetics of a tumour. Detailed understanding of tumour growth kinetics is especially important in co-growing populations, as we here show that the convergence towards an equilibrium—which sets the time scale for potential treatment and relapse effects—sensitively depends on the microscopic cellular growth rates. The often performed, and mathematically convenient rescaling of time that leads to replicator equations might eliminate effects that are crucial for understanding transitions between equilibria and describing relevant time scales of tumour evolution.

## Data accessibility

All data used in this paper are presented in the appendix and can be found online with the cited references.

## Authors' contributions

All authors conceived the study, analysed the data, performed mathematical and statistical modelling, and wrote the paper.

## Competing interests

The authors declare no conflict of interest.

## Funding

P.G. gratefully acknowledges support from the Swedish Foundation for Strategic Research (grant no. AM13-0046) and Vetenskapsrådet (grant no. 2014-6095). P.M.A. gratefully acknowledges support from Deutsche Akademie der Naturforscher Leopoldina, grant no. LPDS 2012-12 for part of this work at Harvard University, and generous support from the Moffitt Cancer Center and Research Institute.

## Acknowledgments

We thank Joel S. Brown (Chicago) and Christian Hilbe (Vienna) for valuable comments and suggestions.

## Appendix A. Fixed points and stability

To investigate the stability of the fixed points of (2.3), we denote the right-hand sides by
A 1and calculate the Jacobian at the fixed point (*u**, *s**)
A 2where subscript denotes partial derivative with respect to *u* and *s*.

**A.1. Boundary fixed points**

For the boundary fixed points we find the following:

At (*u**, *s**) = (0, 0), we find that
A 3with eigenvalues *λ*_{1} = *f*_{1}(0) − *f*_{2}(0) and *λ*_{2} = *f*_{2}(0) − *μ* > 0. The last inequality holds because we assumed a positive net growth rate for both cell types for all *u* ∈ [0, 1]. This fixed point is therefore unconditionally unstable.

At (*u**, *s**) = (0, 1), we find that
A 4with eigenvalues *λ*_{1} = *f*_{1}(1) − *f*_{2}(1) and *λ*_{2} = *f*_{1}(1) − *μ* > 0. Again, the inequality holds because we assumed a positive net growth rate for both cell types for all *u* ∈ [0, 1]. This fixed point is therefore unconditionally unstable.

At (*u**, *s**) = (1, *K*(1 − *μ*/*f*_{1}(1)), we find that
with eigenvalues *λ*_{1} = (*μ*/*f*_{1}(1))(*f*_{2}(1) − *f*_{1}(1)) and *λ*_{2} = *μ* − *f*_{1}(1) < 0. This implies that the fixed point is stable iff *f*_{2}(1) < *f*_{1}(1).

At (*u**, *s**) = (0, *K*(1 − *μ*/*f*_{2}(0)), we find that
with eigenvalues *λ*_{1} = (*μ*/*f*_{2}(0))(*f*_{1}(0) − *f*_{2}(0)) and *λ*_{2} = *μ* − *f*_{2}(0) < 0. This implies that the fixed point is stable iff *f*_{1}(0) < *f*_{2}(0).

**A.2. Internal fixed points**

Internal fixed points exist at points where *f*_{1}(*u**) = *f*_{2}(*u**) for 0 < *u** < 1. The corresponding *s*-coordinate is given by solving d*s*/d*t* = 0 in terms of *u* to get *s** = *K*(1 − *μ*/*f*_{1}(*u**)). The Jacobian at such a point is given by
A 5

To say something about the stability of such a point we need to investigate the signs of the eigenvalues of *J*. We do this by looking at the sign of each matrix entry. For now, we assume nothing about the sign of *f*_{1}′(*u**) − *f*′_{2}(*u**) and instead focus on the other factors in each matrix entry.

First, we see that
A 6Further, we have
A 7
A 8
A 9Here, 0 ≤ (1 − *μ*/*f*_{1}(*u**)) < 1 since *f*_{1}(*u*) > *μ* ≥ 0. This implies that
A 10since both terms are negative. Lastly, we see that
A 11since 1 + 2*μ*/*f*_{1}(*u**) > 0.

This implies that we can write the Jacobian as
A 12where *A* > 0, *B* < 0, *C* < 0 and Δ*f* = *f*_{1}′(*u**) − *f*′_{2}(*u**). The eigenvalues of the Jacobian are given by
A 13Now if Δ*f* > 0 then *A*Δ*f* > 0, and . This implies that *λ*_{1} > 0 and *λ*_{2} < 0, and hence the fixed point (*u**, *s**) is unstable.

If on the other hand Δ*f* < 0 then there are three possibilities, either (i) 4*BC*Δ*f* + *A*^{2}Δ*f*^{2} > 0 or (ii) 4*BC*Δ*f* + *A*^{2}Δ*f*^{2} < 0 or (iii) 4*BC*Δ*f* + *A*^{2}Δ*f*^{2} = 0. If (i) holds then which implies that *λ*_{1,2} < 0. If (ii) is the case then is complex and ℜ(*λ*_{1,2}) < 0. Lastly, if (iii) is the case then *λ*_{1,2} = *A*Δ*f*/2 < 0.

This shows that the stability of the stationary point at (*u**, *s**) is fully determined by the sign of Δ*f* = *f*_{1}′(*u**) − *f*′_{2}(*u**). If Δ*f* > 0 the point is unstable and if Δ*f* < 0 then the point is stable.

## Appendix B. Clonal population doubling rates

Here, all rates are given per day; *in vivo* data taken from Marusyk *et al.* [27].

For LoxL3, we used the following population doubling rates (net growth rates):
B 1and the following percentage of S-phase during cell cycle *χ*:
B 2

For IL11, we used the following population doubling rates (net growth rates):
B 3and the following percentage of S-phase during cell cycle *χ*:
B 4

For CCL5, we used the following population doubling rates (net growth rates):
B 5and the following percentage of S-phase during cell cycle *χ*:
B 6

The distributions shown in figure 2 resulted from all possible pairs of these numbers to calculate *α* and *μ*, equations (3.10) and (3.11).

For generation of the *in vitro* distributions we used normally distributed rates (truncated by 0), with a mean death rate of 0.12 d^{−1} (s.d. 0.0672) and a mean birth rate of 1.32 d^{−1} (s.d. 0.048), adapted from Juarez *et al.* [41].

- Received May 9, 2017.
- Accepted August 31, 2017.

- © 2017 The Authors.

Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.