Extinction rates in tumor public goods games

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 tumor 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 timescales, in particular in co-evolutionary 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.


Introduction
The theory of games was devised by von Neumann and 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]. In order 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]. Behavioral 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,11,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 behavior, the precise measurements of payoffs, as observed from individuals' behaviors, is difficult. Milinski et al. determined all but one payoff parameter precisely, in order to observe tit-for-tat strategies in repeated Prisoner's Dilemma games in fish [16]. Kerr et al. showed that E. coli bacteria can be observed to evolve according to rock-paper-scissors type of interactions, if cellular dispersal is minimal As seen by a recent expansion of interesting theoretical considerations focusing on evolutionary games in biology [17] lies in 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.
Tumor cell populations, including cells of the tumor 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 tumor cell populations [21], evolutionary games have become more widely used as a means to theoretically model tumor evolution, especially after tumor 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 tumor cells produces a 'tumor public good' in form of an insulin-like growth factor [23], in form of glycolytic acid and vascular endothelial growth factor [24], or modeling the dynamic equilibrium between lactate respiration and glycolysis in tumor cells [25]. Such non-autonomous effects between tumor 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 and 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 co-evolving tumor 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 tumor 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 We assume that the linear benefit of the public good arises through growth factor diffusion that occurs on a timescale 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 Our analysis in this manuscript is be 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 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 co-evolving heterogeneous cell population.

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.

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 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 co-evolution of the two cell types: defined for x 1 , x 2 ∈ R + . 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 analyze 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.

Analysis
In order to simplify the analysis of the system (2.1) we apply the following change of variables where 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 defined on u ∈ [0, 1] and s ∈ R + . We note that in the case when s is small compared to the carrying capacity K, such that s/K ≈ 0 the system reduces to and 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.

Fixed points
By solving the equations we 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 Here fixed point 1 and 2 are trivial in the sense that they correspond to a system void of cells. Fixed point 3 and 4 correspond to monoclonal populations and are stable if the resident type has a larger growth rate compared to the invading type.
We note that the stability criteria for the non-trivial fixed points at u = 0 and 1, including potential internal fixed points, are identical with those of the two-type replicator equation with payoff functions f 1 and f 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 are 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 ( du dt , ds dt ) at every point. This implies the condition which is known as the manifold equation [31,14]. By substituting du dt and ds dt from (2.3) and letting s = h(u) we obtain the following equation for h(u) This equation is a non-linear 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 where a i (u) are coefficients that depend on u. We insert this ansatz into Eq. (2.6) and equate powers of µ to solve for the a i 's. We do this for i = 0, 1, 2, introducef (u) = u f 1 (u) + (1 − u)f 2 (u), 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 .
Note 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 withh(u) in (2.3), and we get the following expression (to first order in µ): With the unusual pre-factor that is inversely proportional to the total fitness of the population,f (u), 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.

Results and Discussion
It is often argued that pre-factors 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 tumor-public goods game, in which some cells produce a public good at a cost, rendering a benefit to all cells in the population.

Diffusing public goods game
Autocrine production of growth factors is a common feature of cancer cells, and has previously been modeled 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 assumption 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 modeling [35]. Secondly, additional growth factor provision was shown to be rapid and leading to high levels of tumor 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 tumor evolutionary dynamics, aberrant expression of the cooperative signalling molecule was observed on a tumor 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 we 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 Eqs. (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: This system results from Eqs. (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 Eqs. (2.5). The eigenvalues are and hence the free-rider steady state is stable. For the other fixed point (producers dominate) we have making 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 byh Once 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 Thus, 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 In order 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 intial 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 µ →f (u), performed on the logistic system, implying a never-growing population, in which the death rate equals the average birth rate. that provide growth enhancing public goods to the tumor, most notably in experimental work by Marusyk et al. [27]. There, it could be shown that a mixture of certain clones could not explain tumor outgrowth in vivo by simply using superposition of individual clonal birth and death rates. Rather, synergistic tumordriving 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 immuno-therapy [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 ("tumors 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 tumors 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 about 8 hours long, χ serves as a direct estimate for the percent 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

Timescales of in vivo and in vitro cellular expansions
we calculated the mono-clonal birth rates using Thus, given the net doubling rate r = α − µ, it is possible to estimate the death rate with T S fixed to 8 hours. Data for r and χ are given in Appendix B. Since for both r and χ, several independent measurements were performed, 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], see Figure 2 A. In the in vivo tumor growth experiments, exponential growth was observed within the time frame of 50 to 80 days, at growth rates up to two population doublings per day (net growth rate) [27]. However, in most tumors 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 to the birth-death rate ratios observed in cell cultures, where birth rates often exceed death rates by an order of magnitude (α/µ ≈ 10) [42,43,41,44].
As a notable difference to the previous chapter, here we assume both α 1 = α 2 and µ 1 = µ 2 . Thus, instead of (2.3), we now deal with the more general payoff structure and obtain the following ODEs for frequency of producer cells and total size of the system   [27]. Although net tumor growth was high, death and birth rates were similar in all clones considered. In comparison, we also show in vitro cell line rates, estimated by Juarez et al. [41]. We further used the fact that the IL11 cells are growth factor producers. (B) Using median birth and death rates from the distributions in (A), we measured the ε-fixation time numerically determined using Eqs. (3.13) (defined as the time to reach an ε-neighborhood equilibrium value of u, with ε = 0.001, u 0 = 0.5) and compared it to the ε-fixation time of the standard replicator equation (3.8). Note that we used Eqs. (3.13) for this numerical procedure. For IL11 we used α 1 = 0.684/day and µ 1 = 0.596/day. For LoxL3 we used α 1 = 0.617/day and µ 1 = 0.515/day For CCL5 we used α 1 = 1.214/day and µ 1 = 1.031/day. β = 1, with u 0 = 0.5 and s 0 = 0.01/K. Note here that the peak in ε-fixation time marks the shift from u → 1 to u → 0 as the cost increases; this transition can only occur when producers and non-producers have similar birth and death rates. (C) Comparison of ε-fixation times determined numerically using (3.13) to the analytical approximation (3.19), parameters the same as in (B). and we seek to estimate the time it takes to reach a small ε neighborhood of the equilibrium |u(t) − u * | ≤ ε, shown in Figure 2 B. The combinations 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 We 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 which 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) v 0 (t) = 1 We 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 where 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 producers 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. Eq. (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, 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 analyze 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 itself governed by the apoptotic rate of the dominant clone, compare to Eq.(2.8).

Summary and Conclusions
We here have presented calculations that were concerned with the stability and time to reach a neighborhood of equilibrium points in evolutionary game dynamics between two types of tumor cells. We focused on the dynamics of a tumor public good (tumor 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 need 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. ). An observation that cannot be explained with our analytical approximations is that the extinction times of the logistic growth dynamics tend to be closer to the (short) extinction times of the replicator dynamics for smaller death rate µ 1 (µ 2 fixed). However, the utility of our approximations is supported by the observation that overall variability in ε-fixation times is driven by µ 2 , the death rate of the dominant cell type.
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,48,49,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 tumor dynamics [58,59,60,25,61,23,34,24]. However, the assumption that the population is either at constant size may be limiting, as also recently discussed by Li et al. in the context of co-growing and co-evolving bacterial species [62]. 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 tumor dynamics of in vivo tumor growth showed that this parameter regime might indeed be biologically relevant (Figure 2), even when the tumor 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,43,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 describe both co-growth, as well as co-evolution of two tumor 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 stable manifold. This framework allowed us to examined the degree of the resulting variability in ε-fixation times based on previously measured in vivo tumor 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 non-linear size changes over time, especially during treatment [63,64,65,66], and treatment can shift the evolutionary game [67]. Furthermore, selection mechanism 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 modeling 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 tumor. Detailed understanding of tumor 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 re-scaling of time that leads to replicator equations might eliminate effects that are crucial for understanding transitions between equilibria and describing relevant time scales of tumor evolution.

Author contributions
All authors conceived the study, analyzed the data, performed mathematical and statistical modeling, and wrote the paper.

Data Accessibility
All data used in this manuscript are presented in the Appendix and can be found online with the cited references.

Funding Statement
PG gratefully would like to acknowledge support from Swedish Foundation for Strategic Research (grant no. AM13-0046) and Vetenskapsrådet (grant no. 2014-6095). PMA 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.

Conflict of interest
The authors declare no conflict of interest.

A Fixed points and stability
In order to investigate the stability of the fixed points of (2.3) we denote the right hand sides by: and calculate the Jacobian at the fixed point (u , s ) where subscript denotes partial derivative with respect to u and s.

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 ds dt = 0 in terms of u to get s = K(1 − µ/f 1 (u )). The Jacobian at such a point is given by In order 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.
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.

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) The distributions shown in Figure 2 resulted from all possible pairs of these numbers to calculate α and µ, Eqs. (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/day (SD 0.0672) and a mean birth rate of 1.32/day (SD 0.048), adapted from Juarez et al. [41].