Modeling the evolution of collective synchrony

TL;DR

This paper combines the Kuramoto model with evolutionary game theory to explore the balance between group synchrony benefits and individual advantages of cheating, predicting how populations maintain synchrony despite cheaters.

Contribution

It introduces a novel model integrating synchronization dynamics with cheating and policing strategies, providing insights into the evolution of collective synchrony.

Findings

Cheaters can exploit synchrony for individual gain without disrupting the group.

Limited tolerance for cheaters maintains overall synchrony in populations.

Model predicts population compositions that balance synchronization and individual advantage.

Abstract

Group synchrony in the animal kingdom is usually associated with mating. Being in sync is likely advantageous, as it may help in luring the opposite sex. Yet there are also disadvantages -- such as the homogenization of the group -- which make it harder for individuals to stand-out. Here we address this tradeoff, bringing together the Kuramoto model with concepts from evolutionary game theory. We focus on the existence of self-interested ''cheaters,'' which have been extensively studied in a variety of species. In our scenario, cheating individuals take part in the synchronous group display but position themselves (in terms of phase) slightly ahead or behind the pack. This allows them to enjoy both the group benefit of advertisement and the individual benefit of being unique. But a group can only tolerate a limited number of such individuals while still achieving synchrony. We therefore incorporate a from of ''policing'' into our model: if an individual strays too far from the group's synchronous phase, they reveal themselves as dishonest and are punished. Our model offers testable predictions regarding natural population compositions, and will hopefully spur further investigation into not only how, but also why, natural systems synchronize.

Equations

$$\frac{d}{dt}{\theta }_i={\omega }_i+\frac{K}{N}\sum _{j=1}^N\sin ({\theta }_j-{\theta }_i+{\alpha }_i),$$

$$R{e}^{i\varphi }=\frac{1}{N}\sum _{j=1}^N{e}^{i{\theta }_j},$$

$$b_i=g(R)f_i\left({\theta }_i-\varphi \right)=g(R)\left[\text{att}({\theta }_i-\varphi )\text{tol}({\theta }_i-\varphi )\right].$$

$$u_i=\frac{K}{N}\sum _{j=1}^n\sin ({\theta }_j-{\theta }_i+{\alpha }_i)$$

$$c_i=\frac{1}{T}{\int }_0^T|u_i(t){|}^{2}dt.$$

$$p_i=b_i-\beta c_i.$$

$$n_{\text{co}}^{(g+1)}$$

$$n_{\text{ch}}^{(g+1)}$$

$$r_i\propto \frac{e^{p_i}}{{\max }_j{e}^{p_j}},$$

$$\frac{d}{dt}{n}_{\text{co}}$$

$$\frac{d}{dt}{n}_{\text{ch}}$$

$$\frac{d}{dt}{\theta }_{\text{co}}$$

$$\frac{d}{dt}{\theta }_{\text{ch}}$$

$$k^{-1}\frac{d}{dt}\Delta {\theta }_0=n_{\text{co}}\sin (-\Delta {\theta }_0+{\alpha }_{\text{ch}})+n_{\text{ch}}\sin ({\alpha }_{\text{ch}})-n_{\text{ch}}\sin (\Delta {\theta }_0).$$

$$\Delta {\theta }_0^{\ast }={\alpha }_{\text{ch}}.$$

$$n_{\text{ch}}<\frac{1}{1-\cos ({\alpha }_{\text{ch}})},$$

$${\alpha }_i^{(\text{offspring})}={\alpha }_i^{(\text{parent})}+I_i{ϵ}_i,$$

$$I_i=\{\begin{array}{ll}1& \text{if }i\in \text{mutation set}\\ 0& \text{otherwise},\end{array}$$

$$g(R)=R,$$

$$f(\Delta \theta )=1-\cos (\Delta \theta ),$$

$$p_{\text{co}}$$

$$p_{\text{ch}}$$

$$p_{\text{mean}}$$

$$g(R)=R,$$

$$f(\Delta \theta )=\{\begin{array}{ll}\cos (4(\Delta \theta -\frac{\pi }{4}))+1,& |\Delta \theta |<\frac{\pi }{2}\\ 0,& |\Delta \theta |\ge \frac{\pi }{2}\end{array}.$$

$$g(R)=R,$$

$$f(\Delta \theta )=\{\begin{array}{ll}\cos (4(\Delta \theta -\frac{\pi }{4}))+1,& |\Delta \theta |<\frac{\pi }{2}\\ 0,& |\Delta \theta |\ge \frac{\pi }{2}\end{array}.$$

$$N_i^{(g+1)}=(1-\gamma )N_i^{(g)}+k_0{N}_i^{(g)}(p_i-p_{\text{mean}}),$$

$$n_i^{(g+1)}{N}^{(g+1)}=n_i^{(g)}{N}^{(g)}\left[1-\gamma +k_0(p_i-p_{\text{mean}})\right].$$

$$n_i^{(g+1)}=n_i^{(g)}\left[1+\frac{k_0}{1-\gamma }(p_i-p_{\text{mean}})\right].$$

$$n_{\text{co}}^{(g+1)}$$

$$n_{\text{ch}}^{(g+1)}$$

$${\dot{\theta }}_{co}=n_{\text{ch}}\sin ({\theta }_{\text{ch}}-{\theta }_{\text{co}}+{\alpha }_{\text{co}})-n_{\text{ch}}\sin ({\alpha }_{\text{ch}}),$$

$${\dot{\theta }}_{ch}=n_{\text{co}}\sin ({\theta }_{\text{co}}-{\theta }_{\text{ch}}+{\alpha }_{\text{ch}})-n_{\text{co}}\sin ({\alpha }_{\text{co}}).$$

$${\lambda }_1=0,{\lambda }_2=-n_{\text{co}}\cos ({\alpha }_{\text{co}})-n_{\text{ch}}\cos ({\alpha }_{\text{ch}}).$$

Full text

1 Introduction

Synchrony plays a crucial role in the behavior of many animals, often manifesting as coordinated acoustic or visual signals [9]. In many cases, synchrony is thought to emerge due to the evolutionary advantages it provides, as opposed to simply being epiphenomenal. A well-known example is synchronization in mating displays111Collective courtship displays—also known as “lekking”—need not be synchronous[16, 24]., where coordinated flashing or calling can amplify signals, making them more effective for attracting mates over longer distances (the “beacon hypothesis”[7]).

However, synchrony also occurs in contexts where its benefits are less obvious, raising intriguing questions. For instance, certain species of male fireflies exhibit synchronized flashing, which is thought to help attract females from afar [7]. Yet this raises a paradox: if all individuals flash simultaneously, how can any one of them stand out to a potential mate? Similar behavior is observed in fiddler crabs, where males wave their claws in synchrony, especially in the presence of females. Strikingly, some species even exhibit this behavior in the absence of females, suggesting that the function of synchrony may go beyond reproductive signaling. These examples highlight the need to better understand the evolutionary origins and stability of collective synchronization.

Empirically, getting traction on these questions poses multiple problems. Precise measurement of these behaviors in natural ecological contexts is challenging [15, 25]; automatically detecting where and when mating events occurr adds another layer of complication. Without such data, though, it is difficult to gain insight into which males are sexually selected for.

We thus focus our efforts here on theory, with the hope that predictions will later be tested against field data. Mathematically, the Kuramoto model has long served as a foundational framework for studying synchronization in large populations [17, 30]. When one accounts for the individual cost of synchronization, it becomes natural to frame the problem using the language of mean-field games [12, 8, 3, 31]. In this approach, each individual seeks to minimize a cost function, and the system settles into a Nash equilibrium that balances synchronization benefits and individual costs. However, standard mean-field game approaches usually fall short with regard to long-term trade-offs, as they only provide the Nash-equilibrium for the current population, and do not say anything about how the population will evolve. As such, the resulting equilibrium does not necessarily reveal which synchronization strategies are preferred or stable over evolutionary time.

Evolutionary game theory provides another lens for studying this problem [32, 29], allowing us to examine which synchronization strategies will be favored and adopted at equilibrium. This approach emphasizes the role of payoff structures in shaping strategy. Yet, in many cases, the choice of the payoff matrix lacks a principled justification, limiting the explanatory power of the framework.

Here we combine the Kuramoto model with an evolutionary game-theoretic framework to investigate the emergence and stability of synchrony in animal groups. Inspired by dishonest behavior (“cheating”) found in quorum sensing bacteria [27], fiddler crabs [5] or even firefly “femmes fatales" [19], we imagine a scenario where cheaters exist within synchronous groups. Cheating, in this context, would mean taking part in the display, but positioning oneself as an outlier with respect to phase: that is, signaling with some phase advance or delay relative to the group. This would allow the cheaters to, on the one hand, appear as cooperative (they are taking part in generating the collective signal), but to also stand out—potentially giving them an advantage relative to the rest. We ask if, and how, they could play a part in the synchronous dynamics.

Finally, we also introduce the concept of policing (found, e.g., in bacteria[34] and honeybees[11]). Without policing, there may be no incentive to cooperate, so this is a key ingredient to stabilizing synchronous dynamics on an evolutionary time scale. As cheaters can hamper collective benefits, if other cooperating individuals detect such activity, they may “enforce” the rules (we don’t specify here how this may be achieved; it could also be manifested differently in different species). As such, when cheaters become too obvious (e.g., they call anti-phase to the rest of the group) they become susceptible to detection and punishment.

2 Results

2.1 Models

Our models incorporate two key time scales: a “fast” (within-generation) time scale representing the oscillatory behavior of an individual during their lifespan, and a “slow” time scale representing multi-generation evolutionary processes.

For the fast time scale, we assume that each individual’s behavior follows the Sakaguchi-Kuramoto model[26] (a modified Kuramoto model[30]):

[TABLE]

where $\theta_{i}$ is the phase of individual $\omega_{i}$ is its intrinsic frequency (which we take to be universal), $K$ is the coupling strength, and $\alpha_{i}$ represents the individual’s phase lag strategy. For simplicity, we assume that the strategy of each individual is fully characterized by its phase lag parameter $\alpha_{i}$. A phase lag of zero would mean that the individual dynamics are geared towards in-phase relationships with neighbors. A nonzero phase lag pushes the individual to be ahead or behind its neighbors in terms of phase (though system with identical individuals with nonzero phase lag, e.g., could still yield synchrony). Phase lags were considered in prior works [35], but were treated as a proxy for delays; here we reinterpret the meaning, and suggest that this acts as a built-in strategy to aim to be different.

On the slow time scale, we model natural selection as optimizing a net “payoff” to each individual. That payoff is dependent on both group-level and individual level properties, and is expressed as the difference between a benefit and a cost. The benefit $b_{i}$ for individual $i$ is defined as the product of two terms: $g_{i}$, the expected number of females attracted to the group (determined by group-level synchrony), and $f_{i}$, an individual “strategy function,” which quantifies the value of an individual’s flash timing strategy. We assume $g_{i}=g_{i}(R)$ is a positive increasing function of the order parameter $R$ defined as

[TABLE]

where $R\in[0,1]$ quantifies the level of synchrony, and $\phi$ is the average phase. The strategy function $f_{i}=f_{i}(\theta_{i}-\phi)$ is a function of phase difference $\theta_{i}-\phi$ 222This is the phase of oscillator $i$ relative to the mean group phase, and we will later denote it by $\Delta\theta$., and is taken to be a product of an “attraction” and a police “tolerance” function, as illustrated in Fig 1. It captures both how individual timing relative to the group influences mate attraction, and how standing out from the group might be tolerated by policing333Policing enables the group to act against extreme (or obvious) cheaters, reducing their payoffs as a result.:

[TABLE]

On the cost side, letting

[TABLE]

denote the active control or effort444assuming that oscillating ($\omega$) isn’t part of the effort, as they are all oscillating equally; so here we only address the effort to adjust relative to others. exerted by individual $i$, we define the individual cost $c_{i}$ over a time window $[0,T]$ as

[TABLE]

The net payoff $p_{i}$ is defined as the difference between the benefit and the cost, weighted by a relative cost parameter $\beta$:

[TABLE]

2.2 Group membership dynamics with binary phase lag

We begin with the simple scenario consisting of two groups of individuals, cooperators and cheaters, each of which comprises a corresponding fraction of the population $n_{\textrm{co}}$ and $n_{\textrm{ch}}$ (with $n_{\textrm{co}}+n_{\textrm{ch}}=1$). Then the evolutionary dynamics corresponding to the above payoff function Eq. (5) satisfies:

[TABLE]

where the index $g$ indicates the generation number and $k_{1}$ quantifies the strength of the linear dependence between individual payoff and reproduction chances. (See Appendix A for the derivation of Eq. (6))

The population evolution in Eq. (6) represents expected values. It can be implemented with stochastic effects by defining the relative reproductive rate $r_{i}$ to follow the Fermi selection rule [22] as a function of payoff,

[TABLE]

so that each individual in the next generation has a probability of r_{i}\big{/}\sum_{j}r_{j} to be the offspring of individual $i$.

Taking $k_{1}=k\Delta t$, where $k$ now captures the time scale, we set $n^{(g)}\mapsto n(t)$, and $n^{(g+1)}\mapsto n(t+\Delta t)$ and let $\Delta t\to 0$ to obtain the continuous version of Eq. (6):

[TABLE]

where $p_{\textrm{mean}}=n_{\textrm{co}}p_{\textrm{co}}+n_{\textrm{ch}}p_{\textrm{ch}}$.

We numerically explore Eqs. (8) with the interpretation of cooperators as individuals with a phase lag of $\alpha_{\textrm{co}}=0$ and cheaters as individuals with a phase lag $\alpha_{\textrm{ch}}\neq 0$ (but equal for all of them—they are identical). We emphasize that here $\alpha_{\textrm{co}}$ and $\alpha_{\textrm{ch}}$ are fixed: selection is applied only to the abundances of group memberships. Another way to view this is that the $\alpha$ values can change, but these changes are restricted to exactly two specific values (namely, zero and $\alpha_{\textrm{ch}}$).

Fig. 2(a) presents the results of numerical simulations with $N=1000$ oscillators and illustrates the critical threshold beyond which a population of cheaters can no longer sustain synchrony. Fig. 2(b) shows the dynamical relaxation to equilibrium in a typical case. For persistent populations, the proportion of cooperators and cheaters at equilibrium is always $50-50$ (as showed in Fig. 2(a)). This is due to the fact that the minority population loses its advantage as it grows to become half the population.

Simple analytical results can be obtained in this case. Assuming that at equilibrium each subpopulation is synchronized and all oscillators are frequency locked: cooperators share the same phase $\theta_{\textrm{co}}$ and cheaters share the same phase $\theta_{\textrm{ch}}$, with $\Delta\theta_{0}=\theta_{\textrm{ch}}-\theta_{\textrm{co}}$ constant over time. $n_{\textrm{co}}$ and $n_{\textrm{ch}}$ are the fractions of cooperators and cheaters in the population respectively. Substituting these into Eq. (1), we obtain:

[TABLE]

or, expressed fully in the angle difference variable $\Delta\theta_{0}$,

[TABLE]

Setting the time derivative to zero, we find the equilibrium solution

[TABLE]

Linear stability analysis reveals that this solution is stable when

[TABLE]

which matches the numerical results shown in Fig. 2(a). In Appendix B we generalize these results to any binary choice of $\alpha$.

The key message we take away from this model: It is possible for both cheaters and cooperators to coexist in the population, but there is an upper limit on the number of cheaters the population can sustain. That upper limit is $1/(1-\cos(\alpha_{ch}))$.

2.3 Evolutionary dynamics with phase lag distributions

In this section, we introduce a new aspect: mutation. The inclusion of mutations is to make the evolutionary dynamics more realistic now that we relax the binarization of the phase lags—now $\alpha$ may change across generations. Due to this modification, we need to include the notion of policing. Without policing, everyone in the group will eventually cheat to increase their reproductive success. As a result, the group won’t be able to synchronize, which in turn will lead the population to die out.

With policing, we allow the strategy function $f(\Delta\theta)$ to incorporate information beyond mere attractiveness: it also now indicates the optimal balance between maximizing mate attraction and minimizing police intervention (due to, e.g., female choice constraints or actions by competitive males). We note that the payoff primarily depends555Since payoff $p_{i}=b_{i}-\beta c_{i}=f_{i}(\Delta\theta)g(R)-\beta c_{i}$ where $g(R)$ is uniform over the population and $c_{i}$ is nearly uniform over the population. on the choice of the strategy function $f(\Delta\theta)$, which itself depends only on deviation from the mean phase. Other components of the payoff are approximately uniform across individuals at equilibrium due to frequency locking.

We explore the dynamics of the evolution of the population’s $\alpha$ distribution numerically. For each generation, to compute individual payoffs at equilibrium, we numerically integrate the dynamics given by Eq. (1). After transients have decayed, we use the final phases to compute the payoff of each individual. The reproduction probability for each individual is determined by Eq. (7) given the payoff.

To include mutation, a random subset of the new population undergoes mutation on the phase lag $\alpha$ such that:

[TABLE]

where

[TABLE]

and $\epsilon_{i}\sim\mathcal{N}(0,\sigma_{\text{mutate}})$ represents the mutation perturbation.

In Fig. 3, we initialize the population’s phase lag distribution as a normal distribution with mean $\mu$ and standard deviation $\sigma$. We consider four cases, all under the assumption that the attracted individuals have no preference for phases advanced versus delayed relative to the group (symmetric attraction function). The four scenarios differ in how policing is applied through the tolerance function: (1) no policing (without any tolerance function), (2) asymmetric policing (e.g., only individuals lagging behind the group are penalized), (3) symmetric policing (individuals both ahead of and behind the group are penalized), and (4) strong policing (individuals even slightly away from the mean group phase are penalized to such an extent that there is no net benefit to being out of sync, despite mate attraction being minimized when in sync).

Fig. 3(a) shows results for the “no policing” scenario. There individuals are incentivized to cheat as much as possible to maximize their individual benefit. This leads to a breakdown of coordination, and the population ultimately goes extinct due to the resulting low level of synchrony.

Fig. 3(b) shows results for asymmetric policing, and more specifically, for the case where only individuals lagging behind the group are penalized. Here the distribution of phase lags $\alpha$ drifts in the positive direction, with the mean eventually approaching $\pi/2$, when a breakdown of synchrony occurs and the population goes extinct.

In Fig. 3(c), symmetric policing is imposed on both individuals ahead and behind the group. We observe that the population splits into two clusters and stabilizes.

Interestingly, this symmetry in the policing function appears to be necessary for survival of the population. Asymmetry leads eventually to takeover by cheaters, and thus populations that survive, according to this model, must find a way to impose a symmetric cost on either leading or lagging cheaters.

Finally, Fig. 3(d) illustrates the case of strong policing , where the incentive to cheat is entirely eliminated: the cost of deviating from the group outweighs any potential benefit. In this regime, unsurprisingly, the population remains stable and synchronous; the system is very similar to the case with no cheaters.

We note that the long-term survival of the population depends on the interplay between the choice of the policing function and the initial conditions $\mu,\sigma$. If the policing is too weak (i.e., the gap between the peaks in the strategy function is too wide) and/or the initial population is not collaborative enough (i.e., $\mu$ is too far from 0), the population will fail to survive—see Appendix B.1 for exact limits.

Additional comments on the robustness of our results can be found in Appendix C.

3 Discussion and Conclusions

We have developed a novel approach based on the Sakaguchi-Kuramoto model where we swap the usual heterogeneity of natural frequencies for a different source of irregularity: a distribution of phase lags $\alpha$ (usually taken to be a constant, although see some exceptions: [28, 23, 21, 33, 18, 14]). Phase lags can take on different interpretations, including individual delays (e.g., due to computation times internal to each agent) or, as we’ve considered here, purposeful attempts at misalignment with others. Since multiple species are known to synchronize while lekking (e.g., fireflies[6], crickets[10], frogs [1] and crabs [4]—see [9] for a review on the topic), we think our model has the potential for wide application.

Our main result is the discovery that populations of both cooperators and cheaters can coexist even in a highly simplified model, and we analytically derive an upper limit on the fraction of cheaters. This suggests that natural animal collectives might indeed contain a subset of cheaters. In addition, according to our model, symmetric policing—where straying ahead or behind the pack is equally enforced—is the only stable option with the presence of cheaters. Thus if policing is somehow observed in a natural system, we make the prediction that it will be equally applied to both early and late signalers.

The structural instability inherent in this symmetric system (i.e., slight perturbations to a symmetric policing function would lead to extinction) can potentially be resolved by introducing additional stabilization mechanisms. For example, over evolutionary time scales, only populations with symmetric payoff structures may persist, while those lacking this symmetry are driven to extinction.

It is difficult to distinguish cheaters from cooperators based solely on the phase distribution, as cheaters tend to form one synchronized cluster while cooperators form another in the steady state. This may be because both groups incur the same coupling cost in our model, given the universal coupling strength $K$. However, if cheaters and cooperators experience different coupling costs, it may become possible to differentiate them from the phase distribution, as suggested in [13].

One limitation of our study is that we use a global order parameter $R$ to quantify synchrony, which does not account for clusters of phases. For example, when the population splits into two phase clusters, a perfect cooperator may not respond to the global mean phase, but rather to the mean phase of their local cluster.

In addition to cheating and policing, other mechanisms may contribute to the collapse of synchrony in some species while preserving it in others. Selective attention and local avoidance strategies have been observed in frog choruses, which can also disrupt full synchronization [2, 20].

Future work could focus on generalizations of our model. For simplicity we kept the model in an all-to-all setup. Of course the positions of the males in space (whether in 1-, 2- or 3D) could be important. It is easy to imagine that the relative positioning of cheaters, be it within clusters or more uniformly distributed across the group, could influence the overall dynamics.

We conclude with a hope that our model could offer testable predictions for real animal populations. To that end, carefully defining certain individuals as cheaters is key. Ideally, when analyzing pairwise phase differences, some nearby neighbors will appear to have consistent nonzero phase differences666Note that nonzero phase differences could also arise simply due to distance—minimal imperfections between neighbors can accumulate causing a rippling affect such that farther neighbors may appear to have some larger phase differences.. If indeed cheaters are well defined and detected, then a rigorous quantification of their fraction within the group could be assessed, yielding insight into the real-world utility of our model.

4 Methods

4.1 Binary phase lag simulation

We start with the simple case where we only allow the phase lag to take on two given values: $\alpha_{\textrm{co}}$ and $\alpha_{\textrm{ch}}$, for cooperators and cheaters respectively. The population sizes of cooperators and cheaters evolve according to Eq. (8). The simulation results are deterministic given the initial conditions.

We use the following functions to calculate the payoff:

[TABLE]

where $\Delta\theta=\theta-\phi$. These are simple choices consistent with our assumptions that $g(R)$ be a monotonically increasing function of $R$ and $f(\Delta\theta)$ be a $2\pi$ periodic function that penalizes individuals with small phase deviations.

There are two possible equilibrium states of the population. When $n_{\textrm{ch}}<(1-\cos\alpha_{\textrm{ch}})^{-1}$, all cooperators have the same phase $\theta_{co}$ and all cheaters have same phase $\theta_{ch}$ with $\theta_{ch}-\theta_{co}=\alpha_{\textrm{ch}}$ at equilibrium. In this case, $p_{\textrm{co}}$, $p_{\textrm{ch}}$, and $p_{\textrm{mean}}$ can be calculated explicitly as

[TABLE]

with $\exp(i\phi)=n_{\textrm{ch}}\exp(i\alpha_{\textrm{ch}})+n_{\textrm{co}}$ ($\phi$ defined as in Eq. (2)).

Note that we are free to remove the cost term from the payoff here, since the cost is the same for all individuals in this scenario, and thus removing it does not change the value of $(p_{\textrm{co}}-p_{\textrm{mean}})$ or $(p_{\textrm{ch}}-p_{\textrm{mean}})$.

We do not simulate the case when $n_{\textrm{ch}}>(1-\cos\alpha_{\textrm{ch}})^{-1}$ since the order parameter $R$ approaches zero in this regime, indicating that the population cannot maintain synchrony and thus lacks sufficient fitness for reproduction, ultimately leading to extinction.

To numerically determine the critical boundary, we integrate Eq. (1) using various initial population compositions and values of $\alpha_{\textrm{ch}}$. We sweep $\alpha_{\textrm{ch}}$ from [math] to $\pi$. For each value of $\alpha_{\textrm{ch}}$, we conduct a series of simulations with increasing initial fractions of cheaters $n_{\textrm{ch}}(t=0)$. The critical proportion of cheaters $n_{\textrm{ch}}^{\textrm{(crit)}}$ is identified as the smallest cheater fraction that causes the system to become incoherent at equilibrium. This procedure yields a series of $(\alpha_{\textrm{ch}},n_{\textrm{ch}}^{\textrm{(crit)}})$ pairs and, as shown in Fig. 2, they appear to follow the curve $n_{\textrm{ch}}^{\textrm{(crit)}}=(1-\cos\alpha_{\textrm{ch}})^{-1}$, consistent with the analytical prediction.

4.2 Population dynamics simulation with mutation

In these simulations we allow the individual phase lags to mutate across generations. A key practical distinction is that we can obtain a closed-form expressions for equilibrium payoff in the binary case, but we cannot do so in the model with mutation: payoff needs to be evaluated numerically. For each generation, we numerically integrate the dynamics given by Eq. (1) up to $T_{\textrm{final}}=100$ with $N=1000$ individuals. Without loss of generality, we choose $\omega=0$ and $K=1$.

For payoffs with symmetric policing, we choose

[TABLE]

For payoffs with asymmetric (one-side) policing, we choose

[TABLE]

In both cases we set $\beta=1$ to calculate the payoff. The choice of $\beta$ is not critical as cost $c_{i}$ is nearly identical across individuals. The payoff for each individual is calculated with Eq. (5) averaged over the final $10\%$ of the simulation time.

In our simulations, we randomly select $10\%$ of the population to mutate at each generation and model those mutations via additive perturbations of magnitude $\epsilon_{i}\sim\mathcal{N}(0,\sigma)$ as shown in Eq. (12). We use $\sigma=0.05$, since smaller $\sigma$ helps keep the resultant dynamics smooth. The magnitude of these perturbations controls the rate at which the distribution of $\alpha$ drifts across generations.

Appendix A Derivation of two-group equations

We assume that the expected number of offspring for an individual with payoff $p_{i}$, in a population with mean payoff $p_{\textrm{mean}}$, is linearly proportional to that difference $p_{i}-p_{\textrm{mean}}$. So after $g$ generations, if there are $N_{i}^{(g)}$ individuals descended from an initial progenitor $i$, the number of new offspring (in expected value) will be

[TABLE]

where $k_{0}$ is a constant of proportionality for birth and $\gamma$ determines the fraction of individuals who die in a single generation timescale.

Rewriting in terms of the fractional population $n_{i}^{(g)}=N_{i}^{(g)}/N^{(g)}$ where $N^{(g)}=\sum_{i}N_{i}^{(g)}$, we get

[TABLE]

To enforce that $\sum_{i}n_{i}^{(g)}=1$ in all generations (so $n$ represents a fraction as intended), the ratio $N^{(g+1)}/N^{(g)}=1-\gamma$ is required. Then this simplifies to

[TABLE]

We can define $k_{1}=\frac{k_{0}}{1-\gamma}$ as an effective generational growth rate for each subpopulation $i$.

In the simple case where there are only two types of individuals present in the population, this simplifies to two coupled equations. Considering the two types as cooperators and cheaters, the system becomes

[TABLE]

where $n_{\textrm{co}},n_{\textrm{ch}}$ represents the fraction of the population of cooperators and cheaters respectively, and the index $g$ indicates the generation number.

Appendix B Linear stability analysis of the binary phase lag model

Assume at equilibrium, the whole population is frequency locked. The cooperators all have phase $\theta_{\textrm{co}}$, and the cheaters all have phase $\theta_{\textrm{ch}}$, with $\Delta\theta_{0}=\theta_{\textrm{ch}}-\theta_{\textrm{co}}$ constant over time. $n_{\textrm{co}},n_{\textrm{ch}}$ is the percentage of cooperators and cheaters in the population respectively ($n_{\textrm{co}}+n_{\textrm{ch}}=1$). Without loss of generality, we can set $\omega=0$ and $K=1$. Then we have (in a co-rotating frame)

[TABLE]

The two eigenvalues of the Jacobian matrix for this system are:

[TABLE]

Note that $\lambda_{1}$ is the eigenvalue corresponding to the direction $\theta_{ch}-\theta_{co}=\alpha_{\textrm{ch}}-\alpha_{\textrm{co}}$, i.e., perturbations that preserve the phase difference between cheaters and cooperators; the system is neutral to any perturbation in this direction, which simply reflects the freedom of choice of phase coordinate origin.

Setting $\lambda_{2}<0$, we find

[TABLE]

In particular, when $\alpha_{\textrm{co}}=0$ and $\alpha_{\textrm{ch}}=\alpha$, this reduces to the stability condition

[TABLE]

This sets an upper bound on the proportion of cheaters that may be present in a synchronous population.

B.1 Connection between policing and initial cooperativeness

Here we comment on the connection between the strategy function and the initial population distribution for long-term stability of the population.

Suppose $\sigma_{\text{mutate}}$ is small and the gap between the two modes of the strategy function is $\ell$. We observe that the population eventually stabilizes into two distinct subpopulations. The separation between the two modes of the strategy function approximately matches the phase lag difference between the two stabilized subpopulations, as individuals maximizing their benefits. Because the mutation process has zero mean (no directional bias for $\alpha$), the two subpopulations are expected to be equal in size and symmetrically located about the initial $\alpha$ distribution’s mean $\mu$ (without lose of generality we choose $\mu>0$). That is, $\mu_{1}=\mu-\ell/2$ and $\mu_{2}=\mu+\ell/2$. We focus on the case where $\mu_{1}\in(-\pi/2,\pi/2),\mu_{2}\in[\pi/2,3\pi/2]$, which corresponds to parameter regimes in which the population can persist in the long term.

Assuming that the standard deviations of the two subpopulations are small, the stability condition can be approximated using Eq. 17, with $\alpha_{\textrm{co}}=\mu_{1},\alpha_{\textrm{ch}}=\mu_{2},n_{\textrm{ch}}=1/2$. This leads directly to the condition:

[TABLE]

Appendix C Robustness

C.1 Robustness: discrete vs. continous evolution

For the model with binary $\alpha$, we check that results from simulations with discrete generations (Eq. (6)) are consistent with those from equivalent simulations with continuous evolution (Eq. (8)). In Figure 4, as expected, the discrete simulation appears as a noisier version of the continuous one, with the overall picture remaining consistent.

C.2 Robustness: asymmetric bimodal strategy function

We examine numerically whether a symmetric choice of the policing function is essential for the long-term survival of the population. Specifically, in Fig. 5, we consider an asymmetric bimodal strategy function and observe that the subpopulation with negative phase lag eventually dies out due to the accumulated advantage of those with positive phase lag. As a result, the entire population fails to achieve synchronization.

C.3 Robustness: nonzero payoff for perfect cooperators

We examine the robustness of the simulation results when perfect cooperators ($\alpha=0$) receive a nonzero payoff. As shown in Fig. 6, the results exhibit the same qualitative behavior as those in Fig. 3(c) of the main text, where perfect cooperators receive zero payoff.

C.4 Robustness: natural frequency heterogeneity

We examine the robustness of the simulation results when the population has a heterogeneous natural frequency distribution $\omega\sim\mathcal{N}(\mu_{\omega},\sigma_{\omega})$ with nonzero $\sigma_{\omega}$. As shown in Fig. 7 the results appear to exhibit the same qualitative behavior as those in Fig. 3(c) of the main text, where there is zero heterogeneity in natural frequency.