Directed migration shapes cooperation in spatial ecological public goods games

From the microscopic to the macroscopic level, biological life exhibits directed migration in response to environmental conditions. Chemotaxis enables microbes to sense and move towards nutrient-rich regions or to avoid toxic ones. Socio-economic factors drive human populations from rural to urban areas. The effect of collective movement is especially significant when triggered in response to the generation of public goods. Microbial communities can, for instance, alter their environment through the secretion of extracellular substances. Some substances provide antibiotic-resistance, others provide access to nutrients or promote motility. However, in all cases the maintenance of public goods requires costly cooperation and is consequently susceptible to exploitation. The threat of exploitation becomes even more acute with motile individuals because defectors can avoid the consequences of their cheating. Here, we propose a model to investigate the effects of targeted migration and analyze the interplay between social conflicts and migration in ecological public goods. In particular, individuals can locate attractive regions by moving towards higher cooperator densities or avoid unattractive regions by moving away from defectors. Both migration patterns not only shape an individual’s immediate environment but also affects the entire population. For example, defectors hunting cooperators have a homogenizing effect on population densities. This limits the production of the public good and hence inhibits the growth of the population. In contrast, aggregating cooperators promote the spontaneous formation of patterns through heterogeneous density distributions. The positive feedback between cooperator aggregation and public goods production, however, poses analytical and numerical challenges due to its tendency to develop discontinuous distributions. Thus, different modes of directed migration bear the potential to enhance or inhibit the emergence of complex and sometimes dynamic spatial arrangements. Interestingly, whenever patterns emerge, cooperation is promoted, on average, population densities rise, and the risk of extinction is reduced.

The largest real part of the eigenvalues of Eq. (S3.1) is a function of k and l and 11 called the dispersion relation, λ(k, l). Any mode with λ(k, l) > 0 is amplified, while 12 those with λ(k, l) < 0 are dampened. If λ(k, l) < 0 ∀k, l the homogeneous state is stable. 13 The most unstable mode (k * , l * ) is given by max λ(k, l) = λ(k * , l * ) > 0 and is 14 particularly important because it dominates all other unstable modes due to its 15 exponential amplification. The characteristic length scale of the resulting density 16 pattern is inversely proportional to k * for cooperators and to l * for defectors. Numerical 17 observations indicate that these characteristic length scales match closely, k * ≈ l * , due 18 to the dependence of defectors on cooperators for survival (as d > b). Therefore, we 19 restrict the analysis to perturbations with dependent wavelengths, k = l, and abbreviate 20 λ(k) = λ(k, l). Note that λ(0) > 0 iff Q is unstable, i.e. for r < r Hopf , and hence k = 0 is 21 referred to as the temporal mode because it only triggers changes over time. If λ(0) > 0 22 and λ(k * ) > 0 for some k * > 0, temporal and spatial instabilities interfere and can  Table A. Thresholds for pattern formation. Quantitative comparison of analytical predictions from the dispersion relation and from numerical integration for the onset of pattern formation. For activating directed migration A C , R D the threshold marks a lower bound but an upper bound for inhibiting directed migration A D , R C . Numerical thresholds systematically deviate from analytical predictions by requiring higher (lower) directed migration rates for the onset (suppression) of patterns. This difference is expected to decrease when increasing integration times. Parameters as in Fig. 1 but with D D = 0.5 for A C , R D and D D = 0.7 for A D and R C and t = 2000. Initial configuration: homogeneous densities Q perturbed by Gaussian noise with standard deviation 0.01. result in dynamical patterns. For k * → ∞, the emerging patterns have a characteristic 24 length of 0, which implies that no spatial discretization is fine enough to capture them. 25 The asymptotic behaviour of λ(k) for k → ∞ is determined by J S , Eq. (5b). More 26 specifically, det(J S ) > 0 and tr(J S ) < 0 ensure that the largest eigenvalue has negative 27 real part such that λ(k) < 0 for k → ∞. This is the case if the diffusion of cooperators 28 outweighs their aggregation, D C > A C u eq w eq , at the homogeneous density Q.

29
The impact of separate increases in migration rates A C , A D , R C , R D on the 30 dispersion relation is illustrated in S1 Figure The dominant mode k * increases with A C , which reduces the characteristic length of 33 patterns. Conversely, increasing R D decreases k * and increases the characteristic length 34 scale. In contrast, increasing cooperator flight, R C , and hunting by defectors, A D , 35 reduces λ(k) and eliminates unstable modes. The thresholds for pattern formation 36 predicted by the dispersion relation agree well with numerical integration, see Table A 37 in S3 Appendix.
An equivalent requirement is r > r Hopf with r sufficiently close to r Hopf to ensure that 52 Q is a focus and the activator-inhibitor relation of cooperators and defectors is 53 maintained. From the stability of Q follows that the temporal mode is stable, λ(0) < 0. 54 Additionally, we require λ(k) < 0 for k → ∞ and thus that the entries of the Jacobian 55 J S , see Eq. (5b), satisfy s CC < 0, s CD < 0 and s DC > 0, s DD < 0 or, equivalently, 56 D C > A C u eq w eq .

57
Since tr(J I + k 2 J S ) = a CC + a DD + k 2 (s CC + s DD ) < 0 the stability condition for  det(J I + k 2 J S ) < 0 must hold at its minimum, k * , which is given by At k * the stability condition det(J I + k 2 J S ) < 0 simplifies to  Unfortunately, however, the magnitude of this increase sensitively depends on game 75 parameters as well as migration rates and thus largely eludes intuitive interpretations.

76
Nevertheless, for the impact of migration we have 77 det(J S ) = (D C − A C u eq w eq )(R D v eq w eq + D D ) + u eq v eq w 2 eq A D R C , (S3.6) which means that increases in directed and undirected migration rates all increase the 78 threshold with the exception of increases in cooperator aggregation, A C .

79
In summary, our analysis is based on the following three assumptions: 3. det(J I + k 2 J S ) > 0 for some k > 0.

83
Note that condition 3 is essential for unstable stable modes to induce pattern formation, 84 while conditions 1 and 2 are likely more conservative than necessary. More specifically, 85 the formal approach to determine conditions Turing instabilities requires that Q is 86 stable (condition 1). However, numerical results indicate that the criterion maintains 87 July 15, 2019 4/5 relevance for r < r Hopf and was termed diffusion-induced co-existence as opposed to 88 diffusion induced instability for r > r Hopf [3]. Because of the unstable temporal mode 89 emerging patterns may not be static but rather change intermittently or even exhibit