Snowdrift game induces pattern formation in systems of self-propelled particles

Evolutionary games between species are known to lead to intriguing spatiotemporal patterns in systems of diffusing agents. However, the role of interspecies interactions is hardly studied when agents are (self-)propelled, as is the case in many biological systems. Here, we combine aspects from active matter and evolutionary game theory and study a system of two species whose individuals are (self-)propelled and interact through a snowdrift game. We derive hydrodynamic equations for the density and velocity fields of both species from which we identify parameter regimes in which one or both species form macroscopic orientational order as well as regimes of propagating wave patterns. Interestingly, we find simultaneous wave patterns in both species that result from the interplay between alignment and snowdrift interactions—a feedback mechanism that we call game-induced pattern formation. We test these results in agent-based simulations and confirm the different regimes of order and spatiotemporal patterns as well as game-induced pattern formation.

Figure 1

Illustration of self-propelled dynamics: (a) Particles change their orientation $\theta$ by a stochastic angle $\eta$ to ${\theta }^{\prime }=\theta +\eta$ during self-diffusion events. (b) Binary collisions between particles of the same species lead to polar half-angle alignment of their orientations. We allow for additional stochastic noise terms ${\eta }_1$ and ${\eta }_2$ during the collision event.

Figure 2

Illustration of the snowdrift game. (a) Payoff matrix: There is a nonzero payoff only if individuals belonging to different species meet. Then, $A$ receives a payoff $\tau$, while $B$ receives a payoff $\lambda \tau$. The relative fitness $\lambda$ and game strength $\tau$ serve as control parameters. (b) Flow diagram of the replicator equation, Eq. (7), for species $A$ with population density $a$: Unstable fixed points are marked with an open circle, while the stable fixed point at $a_0$ is marked with a filled circle. The arrows indicate the flow toward the stable fixed point.

Figure 3

Solutions of the spatially uniform equations for the velocity fields ${\alpha }_1$ and ${\beta }_1$ [see Eqs. (21)]. (a) Amplitudes of the solutions as functions of the noise strength $\sigma$ with the relative fitness set to $\lambda =1, \lambda =0.6$, and $\lambda =0.2$, as indicated in the graph. Solid lines are stable solutions, dashed lines unstable solutions. If the noise strength is smaller than a threshold value, $\sigma \le {\sigma }_t^A(\sigma \le {\sigma }_t^B)$, then the disordered solution $|{\alpha }_1^{\left(0\right)}|=0(|{\beta }_1^{\left(0\right)}|=0)$ becomes unstable and a new stable solution with $|{\alpha }_1^{\left(0\right)}|>0(|{\beta }_1^{\left(0\right)}|>0)$ appears. (b) Bifurcation (phase) diagram for the spatially uniform solution in control parameter space ($\lambda ,\sigma$). At ${\sigma }_t^A\left(\lambda \right)$, there is a transitions from a disordered state (fixed point $|{\alpha }_1^{\left(0\right)}|=|{\beta }_1^{\left(0\right)}|=0$) to a broken symmetry state with macroscopic polar order in $A$ only ($|{\alpha }_1^{\left(0\right)}|>0, |{\beta }_1^{\left(0\right)}|=0$). $A$ second transition to a state with polar order in both species ($|{\alpha }_1^{\left(0\right)}|>0$ and $|{\beta }_1^{\left(0\right)}|>0$) occurs at ${\sigma }_t^B\left(\lambda \right)\le {\sigma }_t^A\left(\lambda \right)$. The threshold values for the noise strengths coincide if the relative game strength $\lambda =1$.

Figure 4

Stability analysis of the hydrodynamic equations Eqs. (14) and (18). (a) Dispersion relation of the eigenvalue with the largest real part, $s_{\mathbf{q}}$, for different values of the noise strength $\sigma$, and fixed game strength $\tau =0.2$ and relative fitness $\lambda =0.4$. $s_{\mathbf{q}}$ is plotted for wave vector $\mathbf{q}$ pointing in the direction parallel to macroscopic order, $\mathbf{q}\parallel {\alpha }_1,{\beta }_1$ (denoted by ${\mathbf{q}}_{\parallel }$). The spatially uniform state, Eqs. (22), is unstable when the maximum of the dispersion relation, $s_{\text{max}}$, is larger zero. (b) Parameter regime (bands) of $\lambda$ (relative fitness) and $\sigma$ (noise strength) where the spatially uniform polar-ordered states are linearly unstable against spatial perturbations (i.e., $s_{\text{max}}>0$), for different values of the game strength $\tau$ (shades of blue) as indicated in the graph. The two bands are bounded from above by the threshold values ${\sigma }_t^A$ and ${\sigma }_t^B$ indicating the onset of spatially uniform polar order for species $A$ and $B$, respectively (dashed lines). The width of the bands decreases with increasing game strength $\tau$. White regions indicate parameter regimes where the spatially uniform ordered states are stable against spatial perturbations. (c, d) Linearly unstable regions in ($\tau ,\sigma$) space and ($\lambda ,\sigma$) space, respectively. The color code compares the two density field components ($\delta a$ and $\delta b$) of the eigenvector corresponding to the eigenvalue with the largest positive real part: ${\delta }_a:=\delta a/(\delta a+\delta b)$ which ranges from 0 (the largest linear growth rate corresponds to an eigenvector with $\delta a$ component equal to zero and a nonzero $\delta b$ component) to 1 ($\delta b$ component equal to zero and nonzero $\delta a$ component).

Figure 5

Phase diagram in control parameter space ($\lambda ,\sigma$) for game strength $\tau =0.2$. Colored regions indicate parameter regimes where one finds phases that exhibit traveling polar waves, traveling polar wave phase A (TWA) and traveling polar wave phase B (TWB). White regions indicate regimes showing spatially uniform phases. Regions of linear instability are in yellow-brown (depending on the value of the maximum of the largest real part of the eigenvalue of the Jacobian, $s_{\text{max}}$). Bistable regions at the transitions between phases are coloured in blue. Blue dots are calculated numerically and indicate transitions between spatially uniform phases, bistable regions and traveling wave phases. The blue star indicates the vicinity of a small parameter regime, where the transition between the partially polar-ordered uniform phase and the TWA phase becomes supercritical.

Figure 6

Typical wave profiles for the density fields (solid lines), $a$ and $b$, and the velocity fields (dashed lines), $|{\alpha }_1^{\parallel }|$ and $|{\beta }_1^{\parallel }|$, in the traveling wave phases TWA (a) and TWB (b). Control parameters are set to relative fitness $\lambda =0.4$, game strength $\tau =0.2$, and noise strength $\sigma =0.67$ in (a), and $\sigma =0.47$ in (b). The insets show how the density and velocity fields are related to each other: At every point in space, the velocity fields ${\alpha }_1$ and ${\beta }_1$ (dashed lines in the insets) are approximately equal to the values $|{\alpha }_1(\mathbf{r},t)|\approx |{\alpha }_1^{\left(0\right)}(\mathbf{r},t)|=\sqrt{{\mu }_1^A(a,b)/\xi \left(a\right)}$ and $|{\beta }_1(\mathbf{r},t)|\approx |{\beta }_1^{\left(0\right)}(\mathbf{r},t)|=\sqrt{{\mu }_1^B(a,b)/\xi \left(b\right)}$ corresponding to spatially uniform solutions [see Eqs. (22)] for the respective local densities $a(\mathbf{r},t)$ and $b(\mathbf{r},t)$ (green solid lines in the insets). The waves in species $B$ in the TWA phase seem to be an exception to this observation (inset in the bottom left figure).

Figure 7

Bifurcations of the spatial average of the velocity fields, $|{\overline{\alpha }}_1|$ and $|{\overline{\beta }}_1|$, and the spatial variation parameters, ${\eta }_A$ and ${\eta }_B$, as a function of the noise strength $\sigma$ for a fixed game strength $\tau =0.2$ and for different values of the relative fitness $\lambda$: (a) $\lambda =1$, (b) $\lambda =0.2$, and (c) $\lambda =0.5$. Blue and red lines consist of sets of dots that each indicate a stable state calculated numerically. The dashed lines indicate the analytical results for the spatially uniform amplitudes of the velocity fields, Eq. (22); compare with Fig. 3. Regions shaded in blue and yellow indicate parameter regimes, in which the system is bistable or the spatially uniform states are linearly unstable with respect to spatially heterogeneous perturbations, respectively.

Figure 8

Wave profiles for noise strength $\sigma =0.68$, relative fitness $\lambda =0.4$ and different game strengths $\tau$ as indicated in the graphs. The system is in the TWA phase. (a) For $\tau =1$, species $A$ and species $B$ exhibit polar waves in the density fields, $a$ and $b$, with polar-ordered wave peaks (nonzero velocity fields $|{\alpha }_1^{\parallel }|$ and $|{\beta }_1^{\parallel }|$). (b) For $\tau =0.1$, we observe an increase of wave amplitudes in $a$ and $|{\alpha }_1^{\parallel }|$ as well as a decrease in wave amplitudes in $b$ and $|{\beta }_1^{\parallel }|$. This effect is stronger the smaller the value of the game strength $\tau$ is. (c) For $\tau =0$, the polar wave pattern in $B$ has completely disappeared and species $B$ is in a disordered, spatially uniform state.

Figure 9

Regimes of different macroscopic order and spatiotemporal patterns obtained from agent-based simulations for game strength $\tau =0.4$ and two distinct values of the game theory interaction radius $d_{\text{game}}$ as indicated in the graph. Yellow areas mark parameter regimes with traveling polar wave patterns, where TWA and TWB mark the regions where species A and B, respectively, promote wave formation. These two phases merge with increasing values of the relative fitness $\lambda$. Gray shaded areas indicate parameter regimes where the system is in either a disordered or uniformly polar-ordered state, as indicated in the figure legend. For all transitions between phases there are transition regions in which both phases—a uniform phase (either ordered or disordered, depending on the transition) and a polar wave pattern—can dynamically change (light blue regions). Note that, as explained at the end of Sec. 4b, due to demixing of species inside the waves, the relative size of the parameter regime covered by the TWA and merged wave phase widens with increasing $d_{\text{game}}$. The vertical orange lines in the right panel indicate the three slices in parameter space displayed in Fig. 10.

Figure 10

Polar order (solid line) and velocity variation (dashed line) against noise strength $\sigma$ for species A (top row) and species B (bottom row) and a set of different values for the relative fitness $\lambda$ (a-c), indicated in the graph. There are parameter regimes with spatially uniform states (white regions, corresponding to gray shaded areas in Fig. 9), wave patterns (yellow shaded region), and bistable behavior where both wave patterns and uniform (dis-)order (blue shaded region) stochastically alternate. We observe three qualitatively different regimes. For $\lambda =0.23$ (a) there are distinct TWA and TWB phases, separated by an intermediate phase. The TWA phase shows nonzero polar order for species A (${\alpha }_1>0$) and finite velocity variations ($\mathrm{\Delta }{\nu }_A>0$), while there is neither significant polar order in species B (${\beta }_1\approx 0$) nor velocity variations $\mathrm{\Delta }{\nu }_B$. In contrast, the TWB phase is characterized by polar order as well as velocity variation in both species A and B. The case $\lambda =0.28$ (b) is phenomenologically similar to (a), but in the TWA phase there is now also polar order (${\beta }_1>0$) and finite velocity variations ($\mathrm{\Delta }{\nu }_B>0$) for species B. The intermediate phase in which species A is uniformly ordered and species B disordered (center white strip) is much smaller than in (a). For $\lambda =0.31$ (c), the TWA and TWB phases merge to one broad, single traveling wave phase which exhibits polar order and finite velocity variation in both species.

Figure 11

Density and velocity profiles in the TWB phase for noise strength $\sigma =0.32$ and relative fitness $\lambda =0.26$. (a) Snapshot from the agent-based simulation with particle positions of species A and B marked in blue and red, respectively. Inset: blow-up showing species demixing inside the wave. (b) Corresponding wave profile as a function of position $x$ averaged over the $y$ direction. For the analysis we partitioned the system into 256 channels along the $x$ axis, counted the particles per channel and evaluated the channel densities $a, b$, and $\rho$ (upper plot) and polar order ${\alpha }_1, {\beta }_1$ (lower plot). There is a clear correlation between the density wave and polar order; outside the wave the system shows a low-density disordered phase.

Figure 12

Stability analysis of the hydrodynamic equations in ($\tau ,\sigma$) space. The color code measures how fast perturbations in the velocity field of species $A$ ($\delta {\alpha }_1$) grow compared to perturbations in the density of species $B$ ($\delta {\beta }_1$): $\delta {\alpha }_1/(\delta {\alpha }_1+\delta {\beta }_1)$ which ranges from 0 (only perturbations $\delta {\beta }_1$ grow) to 1 (only perturbations $\delta {\alpha }_1$ grow).

Figure 13

Examples of different channel subdivisions. Systems of size $L\times L$ with channels perpendicular to the $x$ axis (left) and diagonal (right) are shown. Due to periodic boundary conditions channels not aligned with either $x$ or $y$ axis may be split in two parts. Channels always have area $L^2/C$.

Figure 14

Snapshots (top row), density (middle row) and polar order (bottom row) of a set of different noise strength $\sigma$ as indicated in the graph; the relative fitness is fixed to $\lambda =0.23$ and the interaction range to $d_{\text{game}}=10$. (Remaining parameters can be found in Table 1.) The profiles for the densities and the polar order are displayed along the direction of polar order as described in Appendix pp4-s2. In the density plot (middle row), the density of species B is adjusted by a factor of inverse relative fitness ${\lambda }^{-1}$ to emphasize the relation between the snowdrift game and local densities of both species; note that at the mean-field level one has $a=b/\lambda$. (a) Disordered phase ($\sigma =0.71$): The system neither shows any spatial patterns (middle panel) nor polar order (bottom panel). (b) TWA phase ($\sigma =0.65$): Species A forms polar ordered waves (direction of order indicated by blue arrow). Species B forms a density wave mirroring species A but does not show polar order. (c) Uniform polar order of species A $(\sigma =0.56$): The system is spatially uniform with no spatial patterns visible. Species A shows uniform polar order, while species B is still disordered. (d) TWB phase ($\sigma =0.40$): Both species A and B show polar waves (direction of order indicated for each species by an arrow in its respective color). While species B only shows order inside the wave, species A still has order in the wake of the wave, evident by the high polar observable in the region behind the wave. (e) Uniform polar order ($\sigma =0.17$): Both species show polar order over the whole system, no wave pattern is visible.

Figure 15

Size dependence of demixing patterns on $d_{\text{game}}$. (a) Snapshot of a system with slow diffusive motion ($v_0/\tau \ll 1$) $[\lambda =0.4, d_{\text{game}}=10, v_0=0.015]$. Demixing into clusters of species B surrounded by species A is clearly visible. (b) Average radius of spots $r_{\text{spot}}$ as a function of $d_{\text{game}}$ for different values of relative fitness $\lambda$. While the dependence is clearly linear, for larger values of relative fitness $\lambda$, the radius $r_{\text{spot}}$ grows faster for increasing game range $d_{\text{game}}$. Inset: independence of the ratio between species $\delta =b/a$ on $d_{\text{game}}.$