Neighborhood size-effects shape growing population dynamics in evolutionary public goods games

An evolutionary game emerges when a subset of individuals incur costs to provide benefits to all individuals. Public goods games (PGG) cover the essence of such dilemmas in which cooperators are prone to exploitation by defectors. We model the population dynamics of a non-linear PGG and consider density-dependence on the global level, while the game occurs within local neighborhoods. At low cooperation, increases in the public good provide increasing returns. At high cooperation, increases provide diminishing returns. This mechanism leads to diverse evolutionarily stable strategies, including monomorphic and polymorphic populations, and neighborhood-size-driven state changes, resulting in hysteresis between equilibria. Stochastic or strategy-dependent variations in neighborhood sizes favor coexistence by destabilizing monomorphic states. We integrate our model with experiments of cancer cell growth and confirm that our framework describes PGG dynamics observed in cellular populations. Our findings advance the understanding of how neighborhood-size effects in PGG shape the dynamics of growing populations.


Supplementary Figures
Parameters: β = 2, σ = 3, α = 1/day, δ = 0.1/day, κ = 0.5/day, K = 1000.   Number of compartments Normalized population size   Number of compartments Normalized population size 2 Fallacy of averages 6 3 Co-evolution of cooperators and defectors in a public goods game  Let C be the number of cooperators and D the number of defectors. Then, the stochastic process which governs the population dynamics is given by, where We note here that r i (y) = r i (x C /(x C + x D )) is unchanged by considering populations and not concentrations. This can be seen by defining x i = i/K, then Supplementary Materials, page 3 Kimmel et al.: Neighborhood size-effects shape growing population dynamics in evolutionary public goods games The transition rates are then given by with i = C, D. We are now able to write down the master equation Unlike the deterministic model, this stochastic model has absorbing states at C = 0 or D = 0. Hence, states which are unstable in the ODE model are not unstable in the stochastic version.

Passage to the deterministic case
For K 1, it is possible to connect the stochastic model to the deterministic model by making the transformation n C = C/K and n D = D/K. The transition rates become Subjecting the master equation (S4) to a system-size expansion expansion (with 1/K 1) [1], we arrive at a Fokker-Planck equation that describes the associated diffusion process [2] To connect with the deterministic model we now assume K → ∞ and so we keep only the first term. Multiplying both sides of Eq. (S6) by n C and integrating over all possible values. The left-hand side is, It is easy to show that the second term will go to zero since ∂n C ∂n D = 0. A similar computation will result in an ODE for x D , and so we obtain which we recognize as a re-scaled version of Eq. (2).

Stochastic simulations
Now that we have established correspondence between the stochastic and deterministic models, we can investigate the impact of stochasticity in this system. We used the Gillespie algorithm to generate trajectories of the stochastic system.
Simulations were carried out using Julia 0.0.7 according to the following scheme: (a) Initialization: Set biological parameters (e.g. α, β, neighborhood size).
(b) Dynamics loop: (1) Determine random composition of neighborhood seen by a cooperator (N C ) or a defector (N D ).
(3) Determine time till next event (e.g. birth or death of either a cooperator or defector).
(4) Determine which event occurs (based on the transition rates).
(c) Output: Dynamics of the total population (t, C(t), D(t)), Final size (C(t final ), D(t final )) and public good benefit observed.
We determine the neighborhood composition, we draw n individuals from the population using a hypergeometric distribution. In the main text, we took the expected neighborhood composition since it is impractical to compute this each iteration step using an ODE model which assumes continuous population sizes. We see qualitatively similar results between the stochastic and deterministic models for large population sizes (Figures 1, 2). However, the stochastic model admits a nonzero probability of extinction of either population in all cases, which is not possible in the deterministic setting. In addition, a region of coexistence is observed when β = 2, σ = 3 which was not predicted in the deterministic model. This discrepancy can best be explained by the fact that the Gillespie algorithm was only run up to 100,000 events. It is likely in this region of parameter space that a long time might be needed to observe the extinction events predicted by the ODE model.

Neighborhood composition
To obtain a group of n individuals drawn from a total population of size C + D with C individuals are cooperators, we use the hypergeometric distribution. The probability of selecting C = k cooperators in the group is given by In our particular scenario, we are more interested in the probability of selecting n−1 individuals to fill the remaining group and determining the size given that particular individual who is observing the neighborhood. The probability that we have k cooperators in a group given the individual observing is a cooperator is then given by In the stochastic simulations, we draw from these to determine the N C , N D at each iteration. In the deterministic version we took the expected value of these distributions. To do that we consider Eq. (S8) by multiplying by k and averaging over all possibilities. We first make note of the following relationship Supplementary Materials, page 5 Kimmel et al.: Neighborhood size-effects shape growing population dynamics in evolutionary public goods games Using the form in Eq. (S9a) we see that a cooperator observes an expected number of cooperators in the neighborhood given byN (S10) A similar relationship using Eq. (S9b) is observed for a defector The perceived expected proportion is then given by dividing by the neighborhood size, For C, D 1 we see that This implies that the deterministic version has a tendency to underestimate the proportion of cooperation in the neighborhood observed by a defector. In contrast, we see that Hence, This implies that the deterministic version has a tendency to overestimate the proportion of cooperation in the neighborhood observed by a cooperator. For large population sizes, the expected neighborhood composition calculated through the hypergeometric distribution approaches that of the binomial distribution.

Fallacy of averages
In the Main Text, we supposed that the cells see their expected neighborhood, and then calculated growth rates that incorporated these expectation values. However, one may also consider the expected nonlinear public good payoff, that is one could calculate averages at a different point. Mathematically, this is equivalent to comparing Taking the expectation of both sides and noting that Assuming that the population and neighborhood size is sufficiently large, we suppose N D ∼ N [y(n−1), y(1 − y)(n − 1)]. We can quantify the effect of switching the order of operations in Eq. (4b) by analyzing the second term in Eq. (S14), (S15) A similar result is obtained for N C . By abuse of notation, let the main text We now summarize the impact by replacing r i in the main text withr i : • In the absence of producers or cooperators, var(N i ) → 0 and sor i ≈ r i .
• The boundary equilibria stability will be unchanged. This can be noted by observing Eq. (S23)-(S24). The eigenvalues are dependent on r i (0) or r i (1). Again, we would expectr i ≈ r i . Additionally, the time will also be unchanged, however the eigenvectors are modified, which will shift the directions trajectories take to approach the boundary equilibria.
• The location of coexistence points will shift. By noting the coexistence condition Eq. (S19), we can see that (S16) • The stability of coexistence points will shift. By looking at Eq. (S34), we see that the stability condition is dependent on the values ofr i at y = y * .
These results are consistent with those observed in figure 4.

Co-evolution of cooperators and defectors in a public goods game
To begin, we rescale the populations x i → Kx i and so Eq (2) becomes We map the dynamical system from producer/defector population to producer frequency and total population. We define Y = x 1 + x 2 and y = x 1 /Y and Eq. (S17) becomes Based on the form given in Eq. (4), we analyze public good functions of the form r C = α C F C − κ and r D = α D F D , where we set F i (0) = 1, but we note a similar analysis can be done with The correspondence between the two forms is given by the relation G i = α i (F i − 1). We call α i the public-good independent (intrinsic) growth rate.

Equilibria points
The system governed by Eq. (S18) admits a minimum of two fixed points: the boundary equilibria. The (y * , Y * ) are given by Internal equilibria can exist provided that two conditions are met. First, both death rates are either both zero or both nonzero. If δ C = δ D = 0, then coexistence is achieved with Y * = 1 and y is undetermined (a line of non-isolated fixed points). If we suppose that δ C , δ D are both positive, then we arrive at the coexistence condition where Y * = 1 − δi ri(y * ) . Finally we note that all physical trajectories (y ∈ [0, 1] and Y ≥ 0) eventually enter the unit box [0, 1] × [0, 1]. To see this we only need to note that if Y > 1, then both terms in Eq. (S18b) are positive and henceẎ < 0.

Generalized model with K-selection competition matrix
Instead of assuming an identical carrying capacities for all subpopulations, one can introduce a more general version of the rescaled model (S17) [3,4], written as Defining y and Y as before, we arrive at where we defined A the competition matrix and u = (y 1 − y) T the population frequency vector (e.g. (A u) i = A i1 y + A i2 (1 − y)). Note that when a = b = c = d = 1, (A u) i = 1 and the model reduces to Eq. (S18). The boundary equilibria remain, but the coexistence points are now given by a modified version of Eq. (S19), where we see it reduces if a = c and b = d (of which our original model is a special case). We analyze the stability of the boundary points: • Defectors win: Assuming that r D (0) > δ D (as before), we see that d > 1/2 is required for stability. In a similar way, the all-C state is stable if Assuming that r C (1) > δ C (as before), we see that a > 1/2 is required for stability.

Linear stability analysis 4.1 Boundary equilibria
We analyzed the linear stability of the defector-only state first and obtained the Jacobian The eigenvalues of this system are λ 1 = δ D −r D (0), λ 2 = δ D r C (0) r D (0) −δ C . This state is stable provided both eigenvalues are negative, which leads to two conditions: (1) r D (0) > δ D , and (2) r D (0)/δ D > r C (0)/δ C . The first condition simply states the reasonable assertion that the intrinsic growth rate of defectors must exceed its respective death rate. The second condition states that the relative growth rate (the ratio of growth rate to death rate) of defectors exceeds that of producers.
The stability of the producer-only state is analyzed via the Jacobian The two eigenvalues are This state is stable provided that both eigenvalues are negative, which leads to two conditions: (1) r C (1) > δ C , and (2) r C (1)/δ C > r D (1)/δ D . The first condition simply states the reasonable assertion that the intrinsic growth rate of producers must exceed its respective death rate. The second condition states that the relative growth rate of producers exceeds that of defectors. Note the similarity in stability criterion between the two states.

Internal equilibria
What if both states are unstable? To answer this, we consider the function Γ(y) = r C (y) − δ C δ D r D (y). Then we see that the "all C" state is unstable provided that Γ(1) < 0 and the "all D" state is unstable provided that Γ(0) > 0. By the intermediate value theorem, there exists a y * ∈ (0, 1) such that Γ(y * ) = 0, which we recognize from Eq. (S19) as the condition for a coexistence point. This shows the perhaps unsurprising fact that unstable boundary points in this system necessitate the existence of internal equilibria. It does not prove that the coexistence point is stable. Indeed, if there are stable closed orbits, a stable internal fixed point is not required. However, if the impossibility Supplementary Materials, page 9 of closed orbits can be deduced, we can conclude that at least one internal fixed point must be stable. We will show later the impossibility of closed orbits in the case of equal death rates (e.g. δ C = δ D = δ).
It is useful to try to determine the maximum possible number of internal equilibria that can exist given a set of frequency-dependent growth rate function functions r C (y), r D (y). We are interested then in the number of possible zeros of Eq. (S19), the y such that Γ(y) = 0. It is sometimes easier to locate extrema rather than zeros of a function. The correspondence between the number of zeros and extrema can be given by the following lemma: Lemma. Let P be the set of all coexistence points and E be the set of all extrema of the function Γ. The maximum number of coexistence points is given by |P | = |E| + 1.
Proof. It is clear that a horizontal line can cross a function near an extrema at most two times. Ordering the extrema E 1 < E 2 < · · · < E n , we note that the interior extrema intersection points are counted twice. Hence, the number of possible coexistence points is given by |P | = 2|E| − (|E| − 1) = |E| + 1 which proves the claim.
With this lemma, the problem of determining the possible number of coexistence points is equivalent to determining the possible number of extremum of Γ. At this point we must specify a form of the public good function. We then investigate the sensitivity of the number of coexistence points to different public good functions. Extrema are located at the zeros of the function g(y) where we have defined We analyzed the existence of internal equilibria for four types of functions with the following properties • The public good always benefits the population (F i (0) > 1).
• More public good never hurts (F i (y) > 0 for all y).
These three properties describe many types of public goods. A notable exception was studied in the explanation of the Warburg effect, where the products of glycolysis are the public good used by cancer cells and at high quantities are actually harmful to the population [5].
First, consider the "almost identical" public good function r C /r D = A = const. Eq. (S25) becomes Since F i (y) ≥ 0, we conclude that g(y) is bounded away from 0 and hence there are no zeroes of g. This implies no extremum and so by the lemma there can be at most one internal equilibria. A similar result holds for the "always better" public good function. It is clear that if F C (y) > γF D (y) or F C (y) < γF D (y), then g will have no zeros and we can conclude at most one internal equilibrium point is possible. These general results show immediately that the linear public good can never have more than one coexistence point and hence saddle-node bifurcations and other interesting phase diagrams are not possible. Note, these results are independent of the listed properties. A general class of models which satisfies the three properties listed above are sigmoidal functions. We consider two general forms: the Fermi-like function considered in the main text, and the Hill function.
The Hill function where the coefficients are all assumed to be nonnegative and h i ∈ Q is the Hill coefficient. Plugging this into Eq. (S25), and setting to 0 leads to g Hill (y) = − γb 2 d 2 h 2 y 2h1+h2 + b 1 d 1 h 1 y h1+2h2 Supplementary Materials, page 10 which we set to zero. There are three cases h 1 > h 2 , h 1 < h 2 and h 1 = h 2 . Let us investigate the case h 1 > h 2 first. Using Descartes' rule of signs (DRoC), we obtain (−, +, sgn(b 1 h 1 − γb 2 h 2 ), +, −). Thus, DRoC shows that we can have at most four positive zeros if the middle term is positive and two if it is negative. The case h 1 < h 2 is analogous and leads to the same conclusion (with the condition on the middle term being reversed). The case h := h 1 = h 2 is special since this combines many of the terms in Eq. (S28): we require We can thus show that one cannot have alternating signs-hence we can only have at most one positive zero. Suppose the signs would alternate. Then we see that Letting b i = 1 + e σ , d D = e σ , d C = e σ+ β n and h i = β n−1 n is the Fermi function used in the main text. Noting that a monotonic transformation conserves the location of extrema, we immediately see that the properties of the Fermi function actually follow as a corollary from the properties of the general Hill function considered. This is because h i are equal and when h 1 = h 2 we showed that there can be at most two internal equilibria.
The general stability of the internal equilibria is much more complex. The Jacobian of an internal equilibria (y * , Y * ) is given by In this case, it is easier to look at the determinant ∆ and trace τ . The values are A necessary condition for linear stability of the internal equilibria is that ∆ > 0, which is true if r C r D > r D r C , or using Eq. (S19) and (S25), g(y * ) < 0. Let δ = δ 1 = δ 2 . The coexistence points occur at r C (y * ) = r D (y * ), with Y * = 1 − δ r C (y * ) . We note here that for physical solutions we require Y * > 0 which implies that r C (y * ) > δ. The Jacobian is The eigenvalues are λ = 0, (1−2Y * )r C (y)−δ. Based on the previous note of Y * it is clear that the second eigenvalue is always negative and so all coexistence points are at least conditionally stable (there exists a trajectory which approaches the point). However, since one eigenvalue is 0, the nonlinear terms are relevant. We want to find conditions that ensure stability of the coexistence points. Defining u = y − y * , v = Y − Y * and inserting these into Eq. (S18) we obtaiṅ If we are close to the equilibrium (u, v 1), we can replace v = Cu where C is related to the entries in the Jacobian. Plugging this into the above and expanding Γ i in small u we obtaiṅ D are the first nonzero terms in the expansion. Neglecting higher order terms we obtaiṅ If k > l then the state is stable if r

Equal death rates -impossibility of closed orbits
By index theory, we require at least one non-saddle fixed point to be in the interior of a closed orbit. Only interior equilibria can satisfy this requirement (the closed orbit cannot leave the unit box in phase space). Furthermore, if the closed orbit surrounds more than one interior equilibria, it must be an odd number 2k + 1 such that there are k + 1 nodes and k saddles, with k ≥ 0. Consider the special case when δ C = δ D = δ, then Eq. (S18) reduces to The coexistence point condition reduces to r C (y * ) = r D (y * ). It is clear It is impossible for a closed orbit to not transverse the curve Φ. However, the line is also invariant (a trajectory cannot leave once it is on it). Hence closed orbits are not possible, since no trajectory can cross through the line y = y * .

Nonlinear model fit
We used equation (S18a) with the assumption Y K ((1 − Y ) → 1 in the rescaled dynamical system) to perform a nonlinear model fit in Wolfram Mathematica 11.2 (NonlinearModelFit and ParametricNDSolve) to obtain maximum Supplementary Materials, page 12 likelihood estimates of β and σ. We set α = 1/day, as this would only alter the overall time scale and not the measured growth rate differences. We also set κ = 0.25/day, as estimated from Fig. 1 in [6], where under highest growth rate concentration the median difference between non-producer growth rate and producer growth rate amounted to 25% of the intrinsic growth rate. This fitting procedure was applied for every value of FBS concentration in the in vitro growth medium, for every integer value of n between 4 and 40, to result in 'distributions' of values of β and σ. These values of β and σ are reported in Figure 3, panels B and C.
6 Fluctuations in n: quasi-spatial model analysis To simulate the effects of N -compartments subject to periodic mixing, we considered the 2N -coupled system of ODEs given by where i = 1, . . . , N and it is understood that x k,i . Each compartments dynamics are governed by the local interplay of producers and defectors within the compartment. The coupling between compartments is through the carrying capacity K.
Simulations were carried out in MATLAB according to the following scheme: (a) Initialization: Set biological parameters (e.g. α, β, initial population size, number of compartments).
Distribute individuals among the compartments according to multinomial distribution.
(2) If exit criterion is met, break out of the loop.
(3) Redistribute cells according to the multinomial distribution. Go back to (1).
(c) Output: Dynamics of the total population (t, x C (t), x D (t)), Final size (x C (t final ), x D (t final )).
The exit criterion was met when x C x D = 0 for some t or if equilibration was reached. This was done by looking back a given amount of time steps and comparing the mean to the current value. After a set number of times they were within a certain tolerance, the loop would end. To solve the ODE, MATLAB's ode45 built-in solver was used. Aside from the compartments, all biological parameters are assumed fixed in the entire system (e.g. there is no compartment-dependent β). However, a new parameter -the selection time T s is necessary. We investigated the impact of T s for fixed compartment size. An interesting effect is observed in Fig. 1. For T s → ∞, mixing never takes place and in this scenario, the defectors go extinct. For finite T s , one observes coexistence. However, as T s → 0 (time scales of mixing and selection are the same), the population appears to be moving towards fixation.

Critical compartment size
It is clear that N is inversely related to the neighborhood size. Therefore, it is not surprising that there should exist some critical number of compartments at which the producers should be favored. To see this, note that N = 1 is the original case where the neighborhood is the whole population. The defectors are the winners as this is the classic tragedy of the commons problem. On the other side, if N → ∞, each compartment will have at most one individual and hence the neighborhood size is effectively 1. In this case, the producers should win. We then expect there to exist some N crit such that for N > N crit , producers are favored.
To approximate this, we use a mean-field approach where we replace the population size in each compartment by its expected value. We analyze the stability of the all producer population y i = 1 for all i. We require the population size to be Y = 1 − δ C r C,i (1) . We replace Y i by its expected value of Y /N . This collapses the N -compartment model to the original model and we recall the eigenvalues of the "all C" state λ 1 = δ C − r C (1) and λ 2 = δ C r D (1) r C (1) − δ D . Assuming the birth rate exceeds the death rate, the first is satisfied. The critical curve is then determined by λ 2 = 0. Solving this for the critical neighborhood size n crit , we obtain .
Now, the neighborhood size is the population size in the compartment and so its average will be n = Y /N . Plugging this into Eq. (S43) yields (after scaling back carrying capacity) with r C (1) = α C 1 + e σ 1 + e σ−β − κ. (S46)