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

Abstract

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.

Fig 1. Typical 1D simulation leading to producer extinction.

A: Snapshots of the system, represented by concentration of producer cells (U), free-riders (V) and growth factor (G) over time, measured in cell-cycle length (time advances top to bottom). The population game is played in 1D, the panels show the concentrations in space. B: Corresponding trajectory of the average number of producers and free-rider densities in their phase space (U = V on the black dashed line). Due to cell motility, the system reaches the slow manifold (orange-dashed line in bottom panel) fast, and spends most of the time traveling along the slow manifold. The slow manifold was calculated numerically from the well-mixed, ODE model. Dimensionless parameters used: γ_U = γ_V = 0.5, a = 0.9, c = 0.5, r = 1, ϵ = 2 × 10⁻³, β = 5, σ = 2.

Fig 2. The effects of initial conditions and dimensionality.

Comparison over dimensionality and random initial condition. A Transition from random to domain wall initial condition. We let the initial condition ${\overrightarrow{u}}_0=p\overrightarrow{W}+(1-p)\overrightarrow{R}$ where $\overrightarrow{W}$ is the domain wall I.C. and $\overrightarrow{R}$ is the random I.C., hence p can be thought of as a spatial correlation measure. Two lengths L = 50, 100 are shown. The simulation was done in 1D and 50 simulations were done per point. The error bars correspond to plus or minus two standard deviations. B Distribution of extinction times by dimension. Dimensionless parameters used γ_U = γ_V = 0.5, a = 0.9, σ = 2, β = 5, c = 0.5, r = 0.9, ϵ = 2 × 10⁻³.

Fig 3. ODE phase diagrams.

The fixed points are labelled by filled (stable) and hollow (unstable) circles. In all cases one observes the slow manifold which connects the fixed points. Each subplot contains trajectories (green = producers win, red = free-riders win). We also show the impact the nonlinearity has on the shape of the slow manifold. Comparing B, D, we see that the nonlinearity deviates the manifold from a straight line. A phase diagram where producers win. Parameter values a = 0.9, β = 5, σ = 2, c = 0.4, r = 1.15. B phase diagram where free-riders win. Parameter values a = 0.9, β = 5, σ = 2, c = 0.4, r = 1. C phase diagram with bi-stability. Parameter values a = 0.9, β = 5, σ = 2, c = 0.4, r = 1.05. The nonzero death rate c has caused the degeneracy of non-isolated fixed points to collapse, leaving behind a slow manifold along which the dynamics travel. D phase diagram where free-riders win. Parameter values a = 0.9, β = 0.5, σ = 2, c = 4, r = 1.

Fig 4. Extinction and coexistence times.

(A-B) The time needed to reach an ε_exit-radius of the stable fixed point as a function of the non-dimensional death rate c using 100 random initial conditions (the standard deviation was smaller than the point size): A: Extinction time of producers. We observe that the linear approximation (Eqs (12) and (14)) to the slow manifold suffers (black, dashed line) as it does not include the impact of β. The colored curves are given by using the explicit calculation of the time in the well-mixed model (Eqs (8)-(10)). Parameters used: γ_U = γ_V = 0.5, a = 0.9, r = 1, ϵ = 2 × 10⁻³. B: Extinction time of free-riders. The linear approximation (Eqs (15) and (16)) performs very well as it includes the strength of the nonlinearity β. Note that β = 5, 50 (green and red, respectively) overlap. Parameters used: γ_U = γ_V = 0.5, a = 0.9, r = 1.2, ϵ = 2 × 10⁻³. C: Time to reach coexistence state/slow manifold. Time needed to within an ε_exit-radius of the coexistence manifold as a function of the initial concentration of free-riders (the initial concentration of producers is U₀ = 0.01) using uniform conditions. Parameters used: γ_U = γ_V = 0.5, a = 0.9, c = 0, r = 1, ϵ = 2 × 10⁻³. All simulations were performed in the 1D system. Due to the re-scaling of the non-dimensional system, all times can be understood in units of the average cell cycle length (1/α).

Fig 5. Extinction times increase linearly with length, and coexistence times increase quadratically with length.

A: Producer extinction time vs. length of the domain. theoretical lines were obtained using Eq (20) (black, blue), and the red and green lines were obtained using the actual time of the slow manifold from the well-mixed model Eqs ((8)–(10)), with the wave front speed in Eq (19). Parameters used: a = 0.9, c = 0.5, σ = 2.0, r = 0.9, γ_U = γ_V = 0.5, ϵ = 2 × 10⁻³. B: free-rider extinction time vs. length of the domain. theoretical lines were obtained using Eq (22). Parameters used: a = 0.9, c = 0.5, σ = 2, β = 0.5, γ_U = γ_V = 0.5, ϵ = 2 × 10⁻³. C: Analysis of the transition from pulled to pushed fronts. This was only observed when the producer invades. As β increases, we shift from the pulled velocity predicted by linear theory. The transition occurs in the colorless zone. Parameters used: a = 0.9, c = 0.5, r = 1.2, σ = 2, γ_U = γ_V = 0.5, ϵ = 2 × 10⁻³. D: Time to reach coexistence state or time to slow manifold with c = 0. Dashed line obtained from Eq. (S43) (S1 Model Analysis). Parameters used: a = 0.9, σ = 2, β = 0.5, 5, 50, γ_U = γ_V = 0.5, ϵ = 2 × 10⁻³. The times to coexistence for different β values deviate by less than 0.1 cell cycles.