Time scales and wave formation in non-linear spatial public goods games

The co-evolutionary dynamics of competing populations can be strongly affected by frequency-dependent selection and spatial population structure. As co-evolving populations grow into a spatial domain, their initial spatial arrangement and their growth rate differences are important factors that determine the long-term outcome. We here model producer and free-rider co-evolution in the context of a diffusive public good (PG) that is produced by the producers at a cost but evokes local concentration-dependent growth benefits to all. The benefit of the PG can be non-linearly dependent on public good concentration. We consider the spatial growth dynamics of producers and free-riders in one, two and three dimensions by modeling producer cell, free-rider cell and public good densities in space, driven by the processes of birth, death and diffusion (cell movement and public good distribution). Typically, one population goes extinct, but the time-scale of this process varies with initial conditions and the growth rate functions. We establish that spatial variation is transient regardless of dimensionality, and that structured initial conditions lead to increasing times to get close to an extinction state, called ε-extinction time. Further, we find that uncorrelated initial spatial structures do not influence this ε-extinction time in comparison to a corresponding well-mixed (non-spatial) system. In order to estimate the ε-extinction time of either free-riders or producers we derive a slow manifold solution. For invading populations, i.e. for populations that are initially highly segregated, we observe a traveling wave, whose speed can be calculated. Our results provide quantitative predictions for the transient spatial dynamics of cooperative traits under pressure of extinction.


Steady states and linear stability analysis
We begin with Eq. (3) from the main text, which describe the dynamics of the densities of public good (PG) producer cells, U , free rider cells, V , and the public good/growth factor itself, G, all of which are functions that depend on space and time. Our choice to neglect a decay rate for the public good µ G is based on the following observations. First, there are molecules that can serve as public goods, which exhibit low decay rates due to binding and unbinding with cell surface proteins, which enhances persistence of these molecules in the long-term. For example, IGF-II (used by Archetti et al., [1]) was found to exist at low concentrations for an extended period of time in the bounded state [2]. As a first approximation, for these molecules, one can assume that the removal of IGF-II from the system is due to absorption/consumption. It is entirely reasonable that the stability condition depends on the location of G 0 , as can be easily verified by analyzing the linear stability analysis for (0, 0, G 0 ). Furthermore, λ(G 0 ) ≥ 1, and we could observe then that our condition is to simply require that all monomorphic states need to be unstable. As the reviewer points out, a finite molecule-lifetime would collapse this condition to the standard fixed point (0,0,0), via the introduction of an additional (small) decay rate of the form −δG. However, this makes very little qualitative change in all other states, and aside from the collapse of the line (0, 0, G 0 ) to (0, 0, 0), no quantitative change would be observed. For example, in the free-rider only state, the location of the fixed point is unchanged at (0, 1 − cr, 0), the only eigenvalue modified would be −(1 − cr)/ → −(1 − cr + δ)/ . For small public good decay rate δ, this effect will be negligible. Additionally, using the fast-time scale implied by the faster diffusion constant of the PG, one could argue that the impact of this eigenvalue is damped incredibly quickly as it is O(1/ ), where = α/δ is the ratio of cellular growth rate to consumption rate. The producer-only state gains another solution with the additional term. However, this turns out to be a small change as this solution is necessarily negative (unphysical). The physical producer-only state solution maintains stability conditions that are based on the same eigenvalues, but with the location of the state slightly shifted due to the additional (small) δ-term.
In non-dimensional form, these partial differential equations reaḋ The boundary equilibria can be read off by setting either U or V to zero. We obtain . This state is stable if λ(G 0 ) < min(cr, 1 − a + c).
We know that at most one coexistence point can exist and that the point satisfies the relation [3] We can immediately see the unsurprising fact that for coexistence, r > 1. That is, we require µ U < µ V . In our case, we can explicitly solve for G by inserting our form for the nonlinear public good and we obtain Finally, we can use this to solve for U, V by noting that G = U/(U + V ) in the coexistence state. Hence we obtain To analyze each steady state, we need the Jacobian of the spatially-homogeneous version of Eq. (S1)-(S3).
If free-riders win the game, we obtain the eigenvalues cr − 1, c(ar − 1), cr−1 . This state is stable provided that 1 r > max(a, c). If producers win the game, we obtain the eigenvalues 1 − a + c − λ(1), c . This state is stable provided that λ(1) − 1 + a > max λ (1) r , c . The stability of the coexistence state involves a very complicated expression relating all the biological parameters. To gain more insight, we look at the reduced Jacobian obtained by setting G * = U * /(U * + V * ). Analyzing the trace and determinant of the resulting 2x2 matrix shows that in the limit as → 0, the coexistence state is always unstable.
A special case is given by c = 0 with U * + V * = 1 and G * = U * . The Jacobian simplifies considerably and we obtain eigenvalues − 1 , 0, G * (1 − a) − λ(G * ). Hence the coexistence state is stable if G * (1 − a) < λ(G * ). This case is clearly degenerate in the sense that we obtain a non-isolated set of fixed points. It is clear that lim G * →0 λ(G * ) − G * (1 − a) → 1 and so there exists at least some portion of the non-isolated set which is stable when c = 0.
2 Calculating ε-extinction times Consider a general systemẋ = A(x) where A(x) can be a nonlinear function of x and A, x ∈ R N with fixed point x * (e.g. A(x * ) = 0). We expand near the fixed point εz = x − x * with ε 1 (note this is just a general small parameter, and is not the same as = α/δ). This leads tȯ At first order, we haveż = N n=1 ∂A(x * ) ∂xn z n which is simply the linearized equationż = Jz where the components of the Jacobian were given. Assume the point x * is linearly stable and that there exists a largest eigenvalue (closest to 0) λ (i) > λ (j) for all j = i with corresponding eigenvector u (i) . Then as t → ∞, we expect z(t) → cu (i) e λ (i) t . To calculate an approximation to the slow manifold, we can without loss of generality write the eigenvector as (1, a 2 , a 3 , . . . , a N ) where a i ∈ C. This eigenvector approximates the direction of the slow manifold. This leads tȯ One can improve the approximation by including higher order approximations of the slow manifold. However, for our purposes, we only use the linear approximation. We also note that any higher order corrections will only enter for order p ≥ 2. Hence an approximation to the dynamics on a 1D slow manifold is given bẏ For simplicity, we suppose that all terms n ≥ 3 are subdominant and can be neglected. The higher order ones are more complicated, but doable if necessary. We assume the initial value U 0 and we can approach the fixed point U * = 0 (recall that we shifted so that the fixed point is at the origin) exit criterion radius ε exit from either side.
therefore the time to enter the ε exit -radius of the fixed point is given by If a 2 = 0 we have while if a 1 = 0 L'Hopital's rule gives We note that the case a 1 = 0 is similar in nature to a center manifold reduction. Furthermore, some general observations, the final time is inversely proportional to the lowest order (nonzero) coefficient and diverges logarithmically in the slow manifold case, and harmonically in the center manifold case. These relations allow us to read off the time to an event in parameter regions where this approximation is suitable. Consider the approach to the fixed point U * = 0 and suppose a n is the first nonzero coefficient of Eq. (S7), then we see that This provides a useful tool to check the order of the slow manifold. By varying the size of ε exit , we can numerically determine the scaling law to find the first dominant term in the slow manifold.

Slow manifold calculation
We calculate an approximation to the slow manifold explicitly in the case where free-riders win and |c(ar − 1)| 1. We calculate a Taylor expansion of the slow manifold in terms of U , where we assume V = V (U ) and G = G(U ) and we also expand in 1. At first order, we obtain The expression to second order is more complicated, and with the aid of Mathematica, we obtain to second ordeṙ (S14) To leading order, the ε-extinction time is then given by For more details and the calculations for the other parameter regimes and fixed points, see the Mathematica notebook in the "spatialPGG" repository at https://github.com/MathOnco/spatialPGG.

Calculating time to slow manifold/coexistence times
When cell death is negligible, µ U = µ V = 0 (thus c = 0 in the non-dimensional case), coexistence is possible. To see this with c = 0 we can divide Eq. (S1) by Eq. (S2) and arrive at If we make the assumption that G equilibrates rapidly, we make the approximation that (the initial value of the public good). This allows us to solve this equation for U and arrive at with ξ = a−1 λ(G0) . Inserting Eq. (S16) into Eq. (??), we havė We also make the observation that with the range for a ∈ [0, 0.25] we see that |ξ| ≤ 0.25 and so we consider it as a small parameter. For t → ∞, a steady state is given by the solution to Expanding V in ξ 1 we obtain (S18)

Calculation of wavefront speed in 1D
Although traveling wave solutions can not exist indefinitely, they were still observed in our simulations of the invasion processes (free-riders and producers separated by a domain wall) as the primary mechanism to reach a homogeneous state.
Following the work in [4], we seek an approximation to the speed of the wavefront. The instability of the unstable state drives the dynamics. The procedure is as follows. First, we linearize about unstable state, the state that is invaded. We assume an Ansatz, that is of the form of a traveling wave e Λt−Kx = e −K(x−Λ/Kt) . Here it is clear that 1/K is the length scale that determines the width of the transition zone and η = Λ(K)/K is the speed of the wave. The wave travels at the minimum possible speed which is obtained by looking at the extremum dη dK = 0.

Free-rider invasion
Suppose the producer-only steady state is unstable. Then we consider a perturbation about that state, U = 1 − c λ(1)+a−1 + e Λt−KxŨ , V = e Λt−KxṼ and G = 1 + e Λt−KxG , where the tilde variables are meant to be small perturbations. Plugging this into Eq. (S1)-(S3) and neglecting higher order terms in the tilde vaiables, we arrive at a linear system A(Λ, K)p = 0 where p = (ŨṼG) T . We are interested in nontrivial solutions to this linear equation, which implies we are looking for when the determinant is 0. This leads to three separate values for Λ, Since the producer-only state is unstable, it must be that one of the eigenvalues is positive and the relevant eigenvalue is Λ 2 , Since we are interested in the speed of the front, we are looking at Minimizing over K, we arrive at (S19)
In an analogous manner to the producer state, we arrive at Since the free-rider only state is unstable, it must be that ar − 1 > 0 and so the relevant eigenvalue is Λ 1 , Minimizing over K, we arrive at An important side note for the invasion of producers into the domain is in the case that 1 − cr < 0. In this case, the free-rider only state is not only unstable, but is biologically infeasible and so the wave speed calculation will be inaccurate since the trajectory in space is no longer traveling from that point to the producer-only state. Instead, it is traveling from the mass-extinction state towards the producer state. The corresponding eigenvalues Λ can be found in an analogous manner and we obtain the estimate for the traveling wave speed, (S21)

Wavefront speed in higher dimensions
The wave front speed in higher dimensions can be calculated provided that the wave front is strictly 1D. That is, if the wave is of the form u(x, y, z, t) = u(x − ηt), then the wave speed can be calculated and is analogous to the previous section (Table A). In principle, other wavefronts can be calculated, for example, the radially symmetric wave that starts from the center of the 3D domain, which invade outward. However, since the domain is finite, the boundaries play a significant role in affecting the speed and so we restrict our attention to 1D-traveling waves.
Length ε-extinction time (1D, 2D, 3D)  6 Calculating ε-extinction times using a traveling wave ansatz Consider the general class of Fisher models where D, R ∈ R n×n are diagonal matrices containing the diffusion coefficients and growth rates respectively. We introduce two characteristic scales, to be chosen, x/L x = X and τ = t/T , which leads to Define r = max R and D = max D. We will manipulate these different scales to analyze the system. First, suppose This implies that spatial variations can be neglected at large length scales. To balance the remaining two terms, we set T = 1/r and see that This observation implies that the system evolves independently over a large spatial scale and is governed solely by the growth term. Hence after rescaling time back we see that U t = R F ( U ) and so the time for equilibration at this scale is given by τ ODE ( U 0 ), where U 0 is the initial data. Suppose instead that L x D/r. To balance we set T = L 2 x /D. and require that T r. Hence in this regime, we see which implies that at small spatial scales, variations are damped out exponentially quickly on the order of T 1/r. Moving on to the last case, we suppose that L x = D/r. Setting T = 1/r leads to all three terms balancing. At this point, we look for traveling wave solutions of the form U (X, τ ) = U (X − ητ ) = U (z). Inserting this into equation (S23) leads to We observe that the speed of this wave if it exists should have the following characteristic scale η ∼ L x /T = √ rD. Then the time it would take for this wave to travel a distance d = L/2 (half the domain) is given by In our case though, the domain is finite and so this wave will eventually reach the boundary. Will the wave just pass through? No. It is clear that in the case of a no-flux boundary, we require u x (L, t) = U (L − ηt) = 0 for all t, which clearly cannot be satisfied once the wave has reached the boundary. What then is the mechanism to finish equilibration once the traveling wave reaches the boundary? To answer this, we note that it is clear that τ should be a continuous function of L. Furthermore, provided that the wave speed is finite, we expect τ wave propagation (L) to be a linear function of L. Recall also that for L small, we are smaller than the transition width of the traveling wave, and are in the region where spatial variations are smoothed rapidly T 1/r, hence we can replace U by its spatially averaged value u. To see this, consider the average value by integrating over the domain The last part follows where we have defined U = u + L φ and the fact that spatial variations are removed quickly. Thus, as L → 0 it is clear that lim By continuity then It is clear that the last term should tend to 0 if we started with an initial condition which was close to the form of the wave. Despite using a step function initial condition, this time to formation was typically much smaller than the first two terms. We also note that the third term is independent of length provided the domain is sufficiently large enough for wave formation. This allows us to easily calculate numerically the traveling wave speed (see main text for details).
We have shown that we may approximate the total time by this superposition of terms and that the initial data for the ODE term is given by the spatially averaged value of the wave profile, but what is this value? To answer this, we first assume that there is one state U * that is invading the state V * . First, we integrate equation (S22) by transferring into the co-moving reference frame z = x − ηt, the equation for the traveling wave profile from −ξ to ξ in such a way that U (−ξ) = U * , U (ξ) = V * . We will also define for notational simplicity the difference w := U * − V * and the average value µ = (U * + V * )/2. We note that we are actually integrating element-wise and so all vectors appearing in multiplication (or division) are meant to be done element-wise.
We will make use of the following identities Next, we multiply by U and integrate equation (S27) to obtain At this point, we make the observation that equation (S32) is a probability distribution. Hence equation (S35) can be formalized as Where · signifies the average value. We can further refine this by considering once more equation (S31) by multiplying by U and integrating over space to obtain Inserting this expression into equation (S36), we obtain Thus for D/η 1, we may approximate the mean of U as simply the average value between the two states.

Calculation of time to coexistence in the absence of cell death
In the case c = 0 (no death), the ODE model gives an approximate time for coexistence to the line U + V = 1.
Simulations of the spatial model show a similar rapid approach to the line with the average population. However, the system continues to run till it is homogeneous. Since there is no traveling wave in this scenario, the only method to reach homogeneity is that of diffusion. In a similar way to the wavefront problem, we assume that the total time is given by We use the approximation that U + V = 1 for calculating the diffusion time. This uncouples U and V and reduces Eq. (S1)-(S3) toU We expand the solutions of U, V, G using a cosine series U, V, G ∼ a n (t) cos(ω n x), where ω n = nπ L . The solution for U is given by and U 0 (x) := U (x, 0). The solution for V is analogous and we note that the time for G to equilibrate is dependent on U and the time is exponentially small.
As an example, we choose the domain wall initial conditions U = I u H(L/2 − x) and V = I v H(x − L/2). We obtain The exit criterion is satisfied when U t 1 < ε exit and V t 1 < ε exit . Noting that n = 1 governs the time to homogeneity (n > 1 is subdominant), we obtain A similar result is obtained in higher dimensions. The procedure is the same for any initial condition. However, we note that the choice of the exit criterion will change the form of τ diffusion . Examples include U t ∞ or U x p (the p-norm). For example, using U x 1 we obtaiñ It is easy to see that this choice leads to a larger increase for increasing length.

Ruling out pattern formations
We have solved Equations. (S1)-(S3) with different initial conditions in one-,two-and three-dimensional systems.
Though the time to ε-extinction or coexistence can vary, all systems eventually reached spatial homogeneity. A main cause of spatial pattern formations in diffusive systems is through diffusion-driven instability (DDI) which occurs through a mechanism known as a Turing bifurcation. We show here that the model does not admit a DDI via a Turing bifurcation.
To begin, we recall that for DDI to occur, we require the homogeneous fixed state be linearly stable in the absence of diffusion. Stability in the absence of diffusion implies 1 r > max(a, c). For DDI, we require the existence of a perturbation which drives at least one eigenvalue positive. This leads to the condition that either −[c(1 − ar) + k 2 γ U ] > 0 or −(1 − cr + k 2 γ V ) > 0. Clearly, neither condition can be satisfied and so no DDI is possible off the free-rider only state. The situation is analogous for the producer-only state and the extinction state, where the addition of diffusion only further drives the stability of the point and does not cause any instability.
The reasons for the markedly different behavior compared to that observed in [5,6] is primarily due to their system mimicking the behavior of an activator-inhibitor system [7,8], where cooperators play the role of activators and defectors -the role of inhibitors. These systems are well-known to exhibit DDIs. In contrast, our system is based on a competitive Lotka-Volterra model with a growth rate dependent on the amount of public good available, which is determined by the amount of producers.

Different initial conditions converge to the homogeneous solution
Our simulations show that systems with higher initial spatial order take longer to converge to the homogeneous solution, see Figure A. We also observe that even though the total distance the free-rider population must is invade is the same (L/2), the traveling waves in figure A2 occur simultaneously and shows that each wave is traveling L/4. Hence the wave front time should be roughly half. The well-mixed model takes 476.79 cell cycles to fixate. The difference between the two times gives the total time for the two-wave travel, which is 503.51. The speed of the wave is given by Eq. (S19) and is equal to 0.1126. Hence the two waves traveled approximately 56.69, a slight overestimate as it should be 50. The discrepancy is likely due to the interaction of the left and right traveling wave fronts being in close proximity at the beginning. In Figure A3, we show the difference between the mean population trajectory of the traveling wave (both waves travel the same mean path) as compared to the well-mixed (ODE) model, which is the orange, dashed line. As expected the trajectories are not the same in this case (compare with Fig. 1, Main Text) as this is due to the traveling wave following the heteroclinic connection from the stable state into the unstable state [9]. However, note that as we approach the fixed point, we converge to the same path the well-mixed model is taking. This is in agreement with the argument presented in Appendix 6. The trajectory in mean producer and free-rider space. The orange, dashed curve is the slow manifold predicted by the non-spatial model. As expected the spatial average travels a different manifold until it gets close to the fixed point. Other (non-dimensional) parameters used in all panels: a = 0.9, β = 5, σ = 2 r = 1, c = 0.5, γ U = γ V = 0.5, = 2 × 10 −3 .

Uncorrelated initial conditions remove domain-length dependence
To show that uncorrelated initial conditions lead to non-length dependent ε-extinction times, we simulated L = 50, 100, 200, 400 and computed the Moran I statistic with weighting matrix w = δ i,i−1 + δ i,i+1 . The statistic ranges from -1 (perfectly dispersed) to 1 (perfectly separated, and 0 (random). Let y be some spatial data andȳ be the mean, then the Moran statistic is given by where N is the number of grid points and W = i j w ij = 2(N − 1). If there is no spatial correlation, the expected value of this statistic is E  Table B: Moran I statistic measures spatial correlation. Uncorrelated initial conditions as measured by the Moran I statistic, which shows the length of the domain does not impact the time to -extinction when uncorrelated ICs are used.