Equilibria and oscillations in cheat–cooperator dynamics

Abstract Cooperative societies can be threatened by cheats, who invest less in cooperation and exploit the contributions of others. The impact of cheats depends on the extent to which they are maintained in the population. However, different empirical studies, across organisms ranging from RNA replicators to bacteria, have shown diverse cheat–cooperator dynamics. These vary from approaching a stable equilibrium to dynamic cyclical oscillations. The reason for this variation remains unclear. Here, we develop a theoretical model to identify the factors that determine whether dynamics should tend toward stable equilibria or cyclical oscillations. Our analyses show that (1) a combination of both periodic population bottlenecks and density-dependent selection on cheating is required to produce cyclical oscillations and (2) the extent of frequency-dependent selection for cheating can influence the amplitude of these oscillations but does not lead to oscillations alone. Furthermore, we show that stochastic group formation (demographic stochasticity) can generate different forms of oscillation, over a longer time scale, across growth cycles. Our results provide experimentally testable hypotheses for the processes underlying cheat–cooperator dynamics.


General 1.Table of all model parameters
2 Frequency and density dependence (scenario 1)

Dynamics of various weightings of density and frequency dependencies
In Fig. 2 of the main text, we showed the population dynamics quickly become static after carrying capacity has been reached, we tested whether the result was just a special case in this section.As shown in Fig. S1, we found the same dynamics across all combinations of weighting coefficients.These results suggest a relatively static dynamic is prevalent.

Analytical stability analysis
We re-write Eq (2) of the main text as where If α is strictly less than 1, then is a unique, positive equilibrium; it is locally asymptotically stable when all eigenvalues of have negative real part (parameters of G and F have been omitted, for simplicity).The Routh-Hurwitz criteria for two-dimensional systems tell us that the eigenvalue condition on J is satisfied when both its trace is negative and its determinant is positive.Direct calculation shows which implies the equilibrium is locally asymptotic stable.A sketch of the phase plane shows not only that the equilibrium is globally stable but also that periodic (oscillating) trajectories do not occur (Fig S).We formalize the line of reasoning, here, using Dulac's criterion.Define H(N co , N ch ) as Nco N ch and consider the vector field, V , defined (strictly inside the first quadrant of the plane) as where, again, we have omitted the parameters of F and G for simplicity.We calculate the divergence of V to be Because F ′ < 0 and G ′ > 0 we can be sure that ∇ • V < 0. Following Dulac's criterion, then, we formally rule out the existence of periodic solutions to Eq (S1), because ∇ • V neither changes sign nor is identically zero.Because trajectories in the planar system, represented by Eq (S1), are bounded (Fig. S2), the Poincaré-Bendixson Theorem tells us that they must either tend to an equilibrium or a periodic orbit thereof.Given that we have ruled out the latter, it must be true that solutions tend to the only equilibrium we have identified inside the first quadrant.
If α = 1, then we can replace the variables N co and N ch with N = N co + N ch and P = Nco N .With these new variables the dynamics become As an aside, we note that the rate of change in the proportion of cooperators is the product of the trait variance (namely, P (1 − P )) and the change in fitness experienced by an individual who switches from cheating to cooperation (namely, F (N ) − h G(P N ) G(P ) F (N )).The sign of N ′ is solely determined by F (N ).Therefore, over time, N converges monotonically to K (the line N co + N ch = K in Fig 2c).The sign of P ′ is also solely determined by F (N ), and so the evolution of P will follow that of N (P increases when N increases and decreases when N decreases).It follows that P will also tend to some equilibrium in a monotonic fashion.Overall, we rule out the existence of periodic solutions when α = 1.
Figure S2: A sketch of the phase plane associated with Eq (S1).Red arrows show the direction followed by trajectories over time.Solid lines show the set of states for which N' co =0 and N' ch =0, respectively.For the former set, trajectories are instantaneously moving in a purely vertical direction; for the latter set, trajectories are instantaneously moving in a purely horizontal direction.Solid lines intersect at the equilibrium (N * co ,N * ch ), and this equilibrium is asymptotically stable.

Parameters in density and frequency dependence function
Our density and frequency dependence functions share the same the basic structure of where x is density (or frequency) of cooperators in the population, h is the coefficient of relative growth advantage of cheats, which is a proxy for the costs of cooperation, a(orb) is the weighting coefficient of the dependence, s is the shape coefficient adjusting the steepness, and t is the threshold, in which the steepest change in response takes place.The effects of each parameter is shown in Fig. S3: increasing a results in a lower response at low cooperator density (or frequency; Fig. S3a), increasing h results in greater response after threshold is passed (Fig. S3b); increasing s results in steeper response, making cheat's dependence more restricted to a specific range (Fig. S3c); increasing t horizontally shifts the response of cheats (Fig. S3d).

ODE solver checking
We set up the ODE solver (the standard package in Python3: scipy.integrate.odeint) to numerically track the population dynamics.The time in each iteration of the solver is adjustable that bigger leap means less precision but less calculation time.Thus, we checked the precision setting to make sure we are not generating numerical artifacts in our result.As shown in Fig. S4, the dynamics in all panels are close to identical, suggesting all settings have reasonable precision.Thus, we use 20 iterations per time unit for all analysis (Fig. S4b).

Intrinsic growth rate checking (r)
We then checked the effects of intrinsic growth rate to make sure the oscillation is not an artifact from special combinations of growth cycle (i.e., periodic population bottleneck) and growth rate (r).We found the pattern of oscillation is robust as long as r > 1.5 in Fig. S5.Hence, we have been using r = 2 to make sure the results are good representations for most ranges of the intrinsic growth rate.

Intrinsic growth rate, r Relative growth advantage coefficient for cheats, h
Figure S5: Amplitude of oscillation in proportion of cheats across various intrinsic growth rates, and relative growth advantage coefficient of cheat.Brighter colour indicates larger oscillation and the amplitude is calculated from the maximal difference in proportion of cheats in each growth cycle after the dynamics become stabilised.

Effects of the costs of cooperation (h)
We use the coefficient of relative growth advantage of cheats, h, as a proxy for the cost of performing cooperation and we found as h increases, a region of no oscillation expands (arrows in Fig. S6b-d).The region is where both weightings are small, meaning cheats do not depend as much on the density or frequency (proportion) of cooperators.The underlying reason is as h is at a larger value, cheats can out-compete cooperators in a wider range of parameters.
In contrast, the result when weighting is at 1 does not change with h, as cheats cannot eliminate cooperators when there is complete dependence in one of the two weightings.Lastly, we found no oscillation at all when h = 1 as the fitness of cheats never exceeds cooperator in such setting (i.e., cheats are eliminated).

Weighting of density dependence, a Weighting of frequency dependence, b
Figure S6: Amplitude of oscillation across various relative growth advantage coefficient for cheats.Higher h means higher benefit of cheats and brighter colour indicates larger oscillation.

Relative fitness of cheats
One important factor that determines whether cheats or cooperators win is the relative fitness, which is defined as where p ch is the proportion of cheats in the population, w ch is the per capita growth rate of cheats, dN ch N ch dt , and w co is the per capita growth rate of cooperator, dNco Ncodt .When w rel,ch > 1, cheats can outgrow cooperators and potentially eliminate cooperators; when w rel,ch < 1, cooperators grow better than cheats and can potentially eliminate cheats.
We looked into two quantities of relative fitness: (1) the maximal difference within each growth cycle, and (2) the average relative fitness over a single growth cycle.We found the pattern in differences is similar to the pattern in oscillation amplitudes (Fig. S7 and Fig. S6).The only exception is when h = 1 because the brighter colors in Fig. S7a mean cheats are much worse than cooperators, whereas brighter colors other panels mean cheats are doing better than cooperators.This explanation is supported by the average relative fitness (Fig. S8), that we only found relative fitness lower than 1 when h = 1.In addition, we found the region of w rel,ch ≈ 1 expands as h increases because relative fitness is 1 for cheats when there is no cooperator in the population (arrows in Fig. S8b-d).

Effects of dilution ratio (D)
The next factor we investigated is the dilution ratio, which controls the size of population bottleneck between growth cycles.Intuitively, we found when D is too low, there is almost no oscillation (Fig. S9a-b), while the pattern roughly stays the same when D is large (Fig. S9c-d).The reason is when there is insufficient dilution would prevent cooperators from recovering their proportion in the population and eventually lead to extinction (see also Fig. 6a-e in main text).

Effects of growing time (T grow )
Another growth cycle-related factor is the duration of growth time in each growth cycle, which determines whether dilution takes place at lag, exponential, or stationary phase of the growth curve.We found once T grow is greater than 2, the results is identical regardless of the duration of growing time (Fig. S10b-d).This result is also easy to explain as cheats would not have time to grow and increase their proportion if growing time is too short (Fig. S9a).Yet, once cheats have chances to gain proportion, there is little difference in the results because both strains would saturate in densities as growing time increases.

Effects of shape patameter (s d , s f )
Besides the settings in growth cycle, another important aspect of the model is the parameters controlling the density and frequency dependence.We began from the shape parameters, s d and s f , where larger s means steeper response (Fig. S3c).We found increasing s d always help oscillation as yellow area consistently increases as s d increases (Fig. S11).On the contrary, the effect of s f is inconsistent as it helps oscillation when s f < 10 but somewhat hinders oscillation when s f > 10 (Fig. S12).This discrepancy happens because (i) density can fluctuate through dilution processes between growth cycle while frequency does not change through dilution, and (ii) increasing s would narrow the range where response can be different.In other words, the fitness landscape becomes more binary when s increases, while density dependence still have the drive to bounce between lower and higher responses (through dilution), frequency dependence completely relied on relative fitness to oscillate and it is harder to move around when the response is too steep.One final note is we set s d and s f to different values because the range of density is larger than the range of frequency, which is confined in {0, 1}.

Weighting of density dependence, a Weighting of frequency dependence, b
Figure S12: Amplitude of oscillation across various shape parameter of frequency dependence.

Effects of threshold parameter (t d , t f )
The last parameter we examined in scenario 2 is the threshold of the response, t d and t f , which horizontally move the response to different ranges (Fig. S3d).We found changing thresholds result in diverse outcomes.Specifically, increasing t d initially increases the region of oscillation, but when t d > 1.5, increasing threshold hampers oscillations (Fig. S13).This is likely because when threshold is too high or too low, the majority of time in the growth cycle would favour one of the two strains (favouring cheats when threshold is low and cooperators when threshold is high), and prevent oscillation.In one extreme case, when t d = 6 and a > 0.9, there is no oscillation as cooperator density hardly exceeds 6 and cheats cannot grow (Fig. S13d) On the other hand, the pattern of changing thresholds in frequency is similar to changing shapes in frequency dependence (Fig. S14).Following their design, we convert the densities of cheats and cooperators to integer after dilution, and we use binomial process with equal probabilities to re-disperse each unit of cheats and cooperators during repartitioning process.Both processes can help generate stochasticity to the system.

Reducing stochasticity through increasing carrying capacity (K )
In addition to Fig. 4 in main text, an alternative way to manipulate stochasticity is to increase the population size within each subpopulation.This can be done through increasing carrying capacity, K. Nevertheless, the duration of growth cycle and the threshold of dependence have to change because the density of the two strain is going to be at a different level (see figure caption of Fig. S16 for more information).We found the amplitude of oscillation decreases as carrying capacity increases (Fig. S16g-i) and the result is consistent with Fig. 4.

Relative fitness of cheats on a board scale
In Fig. S17, we found the average relative fitness of cheats is above 1 but close to 1 at the regions where oscillating dynamics take place (Fig. 5).Together with Fig. 5, these results suggest oscillation happens at the parameter space (1) where cheats are slightly fitter than cooperators and (2) cheats' fitness is highly variable between subpopulations.

Effects of mixing coefficient (F)
We examined the effects of mixing coefficient as it can potentially change the spatial structure of the entire population and help the spread of cheats or cooperators.We found that greater mixing produces a less periodic dynamic -as the mixing coefficient (F ) increases, the peak signal in harmonic regression become smaller, while the non-peak signals become larger (Fig. S18d& e).This is because faster mixing makes the spread of cheats and cooperators faster, resulting in less periodic dynamics (Fig. S18c& f).4.5 Effects of mixing coefficient when cheats are obligate (F, a or b= 1 ) In addition to Fig. S18, where cheats are facultative that they can grow in the absence of cooperator, we found some interesting result when cheat is obligate.Specifically, we found increasing mixing results in lower proportion of cheats (Fig. S19).Further, increasing mixing produces more frequent oscillations as the peak in harmonic regression moves to the left when F is larger.This is because of the complication that cheats are obligate cheats in this case and cannot grow when cooperator is absent in the group These results suggest mixing is a crucial factor of our model and it can determine the extent of periodicity and noise in the dynamics.

Figure S1 :
Figure S1: Population dynamics across various combinations of weightings in the baseline setting.

Figure S3 :
Figure S3: Parameters of logistic function used in density and frequency dependence terms.
Figure S7: Maximum difference in relative fitness of cheats in each growth cycle.
Figure S9: Amplitude of oscillation across various dilution ratios.
Figure S10: Amplitude of oscillation across various duration of growth phase.
Figure S11: Amplitude of oscillation across various shape parameter of density dependence.
Figure S13: Amplitude of oscillation across various threshold of density dependence.

Figure S15 :
Figure S15: Schematic view of simulation steps in scenario 3, where the simulation codes are modified from Muzuuchi et al, 2022.Following their design, we convert the densities of cheats and cooperators to integer after dilution, and we use binomial process with equal probabilities to re-disperse each unit of cheats and cooperators during repartitioning process.Both processes can help generate stochasticity to the system.

Figure S16 :
Figure S16: Increasing carrying capacity within each subpopulation decreases stochasticity and produce smaller oscillation in proportion of cheats.Growth cycle and threshold are changed to different values to allow population having sufficient growing time (T grow is 8 and t d is 6 for K= 200 ; T grow is 16 and t d is 9 for K= 400 ; a= 1 and b= 0 ).

Figure S17 :
FigureS17: Relative fitness of cheats.Brighter color means cheats having greater fitness over cooperators.Relative fitness is calculated at subpopulation level, based on the proportion of cheats within each subpopulation at the beginning of focal growth cycle.The absolute fitness used for relative fitness is approximated from per capita growth rate, the change in cheat or cooperator density within subpopulation, divided by T grow and the initial densities.We then averaged the relative fitness over all subpopulations and time, and rescaled the relative fitness by dividing it with the relative fitness of cooperators because the raw values are deflated by empty subpopulations.Each grid is the average of 10 repeated simulations.All simulations last 200 growth cycles and data was collected from cycle 51 to 200.The number of subpopulations is set to 4000.

Figure S18 :
Figure S18: Local mixing is crucial for oscillation across growth cycles.(a-c) Population dynamics under default mixing (a, F=1.3 ), twice times of mixing (b, F=2.6 ), and global mixing after the dilution in each growth cycle.In global mixing setting, subpopulations are re-initialised through Poisson sampling after summing the density of cheat and cooperator and after periodic population bottlenecks, so stochasticity is not reduced to zero.(d-f) The corresponding Harmonic regression on proportion of cheats in each setting.Density and frequency dependence are set to the combination where largest coefficient of variation was observed inFig.5 (a=0.9,b=0.25  ).

Figure S19 :
Figure S19: Increasing mixing coefficient results in more frequent oscillation when cheats can only grow when cooperator is present (a=1, b=0 ).

Table S1 :
The list of parameters in the model.Sections column provides the relavent section numbers in this document for each parameters.