Eco-evolutionary Logic of Mutualisms

Abstract

Mutualistic interactions among members of different species are common, seemingly stable, and thus apparently enduring. This is at odds with standard mathematical models based solely on between-species interactions, which show mutualisms to be inherently unstable. Models incorporating parameters for punishment and reward strategies demonstrate that the range of conditions over which stability is observed can be extended; however, the role of community-level dynamics impacted by within-species interactions remains relatively unexplored. Here we develop a general and readily applicable approach for analysing a broad range of mutualisms. By incorporating within-species interactions, we show that mutualisms can be stably maintained across diverse environmental conditions without introducing changes to between-species interaction parameters. Further, a balance of within- and between-species interactions is sufficient to allow the persistence of mutualisms encountering ecological perturbations. Our simple and robust framework resonates with emerging empirical data highlighting the role of community-level interactions and population dynamics in maintaining mutualisms.

1 Introduction

The study of mutualistic relationships—interspecific interactions that benefit both species—is a rich field of theoretical and experimental study [3, 6, 10, 11, 38, 49, 70]. Intuition arising from the short-sighted nature of selection, combined with simple game-theoretic models, calls into question the enduring nature of mutualisms: selection stands to reward species that exploit partner-provided resources, even though exploitation risks the destruction of among-partner interactions [77]. Theoretical models—particularly those that assume extremes of actions, for example, agents that either cooperate, defect, provide or withhold resources—lead to similar conclusions [5, 6].

However, the expectation that mutualisms are short-lived and prone to failure is out of kilter with experimental studies that draw attention to the durability of mutualisms [16, 51, 84]. Important in recent work has been the realisation that mutualisms are sometimes stabilised by interactions with a third species [67]. For example, the mutualism between fungus-growing ants and the fungus they cultivate is stabilised by a third species that parasitises the fungus [55]. Similarly, the mutualistic interaction between figs and fig wasps is stabilised by a third parasitic wasp species that competes directly with pollinators [20]. On occasion, one mutualistic partner is the environment of the other. This is particularly so in the case of interactions between animals and associated microbiomes [47] with the recognition that microbial communities encompass a diversity of interacting types and effects [44, 53]. These examples highlight higher-order interactions within and between species, making studying isolated interactions between species an exceptional case.

Advances in biological understanding have fuelled the development of new theoretical frameworks [24, 62]. Of particular note have been models that change simple game theoretic dynamics by inclusion of higher-order factors such as sanctions and rewards [4, 25, 29]. Unfortunately, these models are often specific to particular systems and lack generality when considering the varied ecologies of the different systems. General models that do exist do take into account interspecific and intraspecific interactions but so far neglected external environmental factors such as phenology [4, 23, 25]. Desirable are general models that take into account the dynamic nature of interactions, and especially models that embrace feedback between ecological factors and ensuing evolutionary responses at both population and community levels [2, 42, 48].

Here we examine the effect of within-species interactions on the maintenance of mutualisms. While a multispecies community factor would be a welcome addition to the analysis via evolutionary games, we focus on two species and devote our attention to realistic intricacies even in this simple system. Using a multi-player evolutionary game theoretic approach [34, 69], we show that within-species interactions exert a powerful effect on the dynamics of between-species interactions and can prevent precariously balanced mutualisms from descending into parasitism. We then use our framework to explore two particular ecological factors in mutualistic scenarios—seasonality of interactions and density-dependence.

Typically, the ever-changing nature of environments constantly challenges the stability of mutualistic relationships. While mutualisms elevate the collective productivity of the involved partners, they are not without risks, as their destinies become intricately intertwined [51]. The dynamics of these interactions remain highly responsive to environmental and ecological shifts, whether stemming from natural processes or human-induced activities [1, 13, 32, 46]. In the case of seasonal or episodic mutualisms, there is a significant potential for phenological partner mismatch due to various environmental factors [72]. Specifically focusing on the domain of mutualism-pollination and the associated phenological mismatch, recent investigations have explored social dilemmas and the impact of seasonal fluctuations of the game parameters within isolated populations [33].

Population sizes are inherently dynamic and subject to fluctuations over time. In scenarios where ecological changes occur rapidly enough to be averaged, their influence on evolutionary dynamics is often disregarded. However, the occurrence of evolution unfolding at remarkably rapid timescales, comparable to those of ecological dynamics, is well documented [8, 37, 71, 79]. Consequently, it becomes imperative to embrace eco-evolutionary dynamics, intertwining ecological and evolutionary processes also for mutualisms [17, 59].

Considering the above two ecological factors, we show that mutualisms can be robust to the destabilising effects of seasonality and density-dependence. Our findings indicate that mutualisms, while often prone to failure at the level of individual interactors, can remain stable due to dynamic feedback between ecological and evolutionary dynamics.

2 Model and Results

2.1 Dynamics of Interspecies Interactions

While indeed a general approach for symbiosis can be adopted [90], we focus specifically on mutualisms. The evolutionary origin of mutualisms is often considered in light of theories concerning the evolution of cooperation supported by standard game theoretic models [2, 19]. We are interested in understanding how the benefits of mutualism are distributed among interacting species [9]. In this context the snowdrift game is a useful metaphor [30, 61] in preference to iterated versions of the prisoners dilemma [81]. In a snowdrift game, partners that cooperate are always tempted to defect. However, it can pay to cooperate even when a partner cheats. Contributions from both species generate a common benefit, but costs may differ among the interacting partners. Thus inclination to cooperate (and the temptation to defect) may have different strengths for the two species. In what follows, the potential for mutualistic interspecies dynamics is proxied via the snowdrift game [9, 35, 83], which is discussed in detail in the Appendix.

Each species is assumed to consist of two types of individual: Generous G and Selfish S (Fig. 1). If most individuals are Generous, then selection will favour Selfish types that either withhold some contribution, or take unfairly of resources from the interacting partner. However, all individuals in the game lose if there is an insufficient number of Generous types, and hence neither species can afford to be completely Selfish. Selection thus favours individuals that act selfishly, while at the same time eliciting generous responses from the other species. The fitness of each of the types within a species depends on the frequency of Generous and Selfish types of the other species. If each individual interacts with only one other individual from the other species then the fitness is a linear function of the frequencies of the two types. However if each individual can interact with a number of individuals from the other species—a more likely scenario [85]—then fitness becomes a non-linear function of the frequencies of the two types. This non-linearity is readily captured by multiplayer rather than two player evolutionary game theory [12, 34, 36]. Multiplayer games have the potential to demonstrate rich dynamics that are qualitatively different than the dynamics generated by linear interspecies games [9, 35].

Fig. 1

Evolutionary dynamics with combined inter-intra-species dynamics. Between-species mutualism is described by the snowdrift game [9, 35, 83]. Species 1 plays a $d_1^{inter}$ player game with species 2 while species 2 plays a $d_2^{inter}$ player game with species 1. Each species has two types of players, Generous and Selfish, that interact both with each other and with members of the other species. For intra-species interactions a general framework of synergy and discounting is applied that generates classical outcomes of evolutionary dynamics [21, 39, 63]

We use 1G and 1S to denote individuals of species 1 that are Generous and Selfish and similarly for the other combinations 2G,  and 2S. The frequency of Generous and Selfish types within in a population are given by $x_{.}$, so for example, $x_{1G}$ is the frequency of 1G. Since it is a frequency of a type within species 1, we have $x_{1G} + x_{1S} = 1$. ${\textbf{x}}_1$ represents the vector $(x_{1G}, x_{1S})$, and similarly for species 2. The fitness of 1G determined by its interspecific interactions can now be denoted $f^{inter}_{1G}({\textbf{x}}_2)$, which is a function of the frequencies of traits in the other species (${\textbf{x}}_2$). The fitness of 2G is thus $f^{inter}_{2G}({\textbf{x}}_1)$, and similarly for the Selfish types. Finally, we use ${\textbf{x}}$ with no subscript to indicate a full description of frequencies across both species: ${\textbf{x}} \equiv ({\textbf{x}}_1, {\textbf{x}}_2)$.

2.2 Dynamics of Intraspecies Interactions

For the ensuing analysis of intraspecies dynamics we use a general multiplayer evolutionary game framework [34, 36]. Figure 1 depicts two interacting species, within which there exist Generous and Selfish individuals (dark and light circles, respectively). For the sake of nomenclature we assume that the trait affecting within-species interactions is the same trait as that affecting interactions between species (although the actual interactions can differ, as we show in Appendix A). Throughout the main text of this article we use a coexistence game to capture intraspecies interactions (shown by $\rightarrow \bullet \leftarrow$). It is possible that the Generous and Selfish types as defined in the interspecies interactions are in fact antagonistic competitors or commensals. Such diverse intraspecies interactions and their resulting dynamics are studied in Appendix A.

Models derived from a population genetic perspective that include demographic dynamics have been proposed before [26, 27] and within species interactions can be interpreted in terms of relatedness [2, 89]. If the two types are Generous and Selfish regarding their interaction at the interspecies level, then what type of intraspecific interactions play out within each species? The cost benefit framework described in [21] and [39], makes possible a transition between the four classic scenarios of evolutionary dynamics [64]: neutrality, dominance, coexistence and bistability. For example, coexistence is possible if both types can invade from rare, and bistability is possible if neither pure strategy can invade the other. For within-species interactions, the fitness of a 1G is denoted $f^{intra}_{1G} ({\textbf{x}}_1)$ and that of 2G is $f^{intra}_{2G} ({\textbf{x}}_2)$, and similarly for the Selfish types.

Evolutionary games typically focus on intraspecies dynamics. Strategies usually evolve within a population and their fate over time is decided using the standard replicator equation [43]. In our study, the fate of strategies is determined jointly by intra as well as interspecies interactions, and the interplay between these two becomes crucial.

2.3 Combined Dynamics

Combining both intra and interspecific dynamics provides a complete picture of all possible interactions. While our focus is solely on mutualisms at the level of the interspecies interactions, within each species there are four possible interactions [39, 64] (dominance of either type, coexistence or bistability). This leads to sixteen different possible combinations, since the character of within-species interactions do not need to be the same for the two different species. Assuming additivity in the fitnesses of inter and intraspecies fitnesses, the combined fitness of each of the two types in the two species is given by

$$\begin{aligned} f_{1G} ({\textbf{x}})= & {} p \, f^{inter}_{1G} ({\textbf{x}}_2) \; + \; (1-p) \, f^{intra}_{1G} ({\textbf{x}}_1) \end{aligned}$$

(1)

$$\begin{aligned} f_{1S} ({\textbf{x}})= & {} p \, f^{inter}_{1S} ({\textbf{x}}_2) \; + \; (1-p) \, f^{intra}_{1S} ({\textbf{x}}_1) \end{aligned}$$

(2)

$$\begin{aligned} f_{2G} ({\textbf{x}})= & {} p \, f^{inter}_{2G} ({\textbf{x}}_1) \; + \; (1-p) \, f^{intra}_{2G} ({\textbf{x}}_2) \end{aligned}$$

(3)

$$\begin{aligned} f_{2S} ({\textbf{x}})= & {} p \, f^{inter}_{2S} ({\textbf{x}}_1) \; + \; (1-p) \, f^{intra}_{2S} ({\textbf{x}}_2) \end{aligned}$$

(4)

The parameter p tunes the relative impact that each of the interactions has on the final fitness. For $p=1$ the well-studied case of Red King dynamics is recovered [35], while for $p=0$ the dynamics of the two species are decoupled and can be individually studied using the synergy/discounting framework of nonlinear social dilemmas [39]. Of interest here is the continuum described by intermediate values of p: we need to track the qualitative dynamics as p changes gradually from 0 to 1, for each of the sixteen combinations (Appendix C). The time evolution of Generous types in both species is then given by the difference between the average fitness of the type as given above and the mean fitnesses of the species ${\bar{f}}_1 = x_{1G} f_{1G} + x_{1S} f_{1S}$ (omitting the functional forms of fitnesses) and similarly for species 2:

$$\begin{aligned} {\dot{x}}_{1G}= & {} r_1 \; x_{1G} \; \left( f_{1G} - {\bar{f}}_1 \right) \end{aligned}$$

(5)

$$\begin{aligned} {\dot{x}}_{2G}= & {} r_2 \; x_{2G} \; \left( f_{2G} - {\bar{f}}_2 \right) . \end{aligned}$$

(6)

Rates of evolution are central to co-evolutionary interactions [78]. In antagonistic interactions such as host-parasite dynamics, it is favourable to evolve faster than the other species so as to persist [80]. However in mutualistic interactions it has been proposed that slower rates of evolution could be favourable [9]. This is so because if two species evolve at different rates (here $r_1$ and $r_2$), then it is possible that the balance of benefits skews in favour of the slower evolving species [9]. Furthermore the balance can also be affected by the number of interacting partners [35]. Variability in the number of players can allow incorporation of a multitude of ecological factors into the analysis of interactions. For example, the number of players involved in a game can be different for each interaction, namely $d^{inter}_1$ and $d^{intra}_1$, and similarly for species 2. Here interspecies interactions are proxied by the multiplayer snowdrift game, which incorporates threshold effects (Appendix A). For example a certain number of Generous cleaner fish may be required to clean their host, or a certain number of Generous ants maybe required to protect lyacenid larva from predators. Since the interaction matrices for inter and intraspecies dynamics are completely different, in principle, it is possible to have different costs and benefit functions for the four games (two snowdrift games from the perspective of each species and the intragames within each species). Here we restrict our attention to the same interspecies game and a variety of intraspecies dynamics (Appendix C).

Interspecies interactions, by themselves, can result in Red King effects and other possible complexities as discussed recently [31]. However inclusion of intraspecies dynamics can generate very different outcomes. Even holding other parameter values the same (such as cost-benefit values, thresholds, number of players etc.) and assuming that intraspecies dynamics account for only one third of cumulative fitness, the qualitative dynamics can be radically different and even result in the stable persistence of Selfish types at intermediate frequencies of types (Fig. 2). For different types of intraspecies interactions, a rich set of possible dynamics emerges (Fig. 6).

2.4 Effect of Seasonality on Interaction Dynamics

Typically, environments fluctuate, providing a further challenge to the stability of mutualisms. While mutualisms raise the combined productivity of their partners, they come with risk, since fates are bound together [51]. Interactions are sensitive to environmental and ecological changes that might occur naturally or be catalysed by anthropogenic activity [1, 13, 32, 46]. A paradigmatic example of such seasonal or otherwise episodic mutualisms are pollination interactions. They run a high risk of phenological partner mismatch as a result of environmental factors [72]. Focusing on social dilemmas, seasonal changes in the game parameters have been recently analysed in single populations [33]. For two populations, this means that the effect of interspecific interaction changes over time. To this end, our model also incorporates seasonality.

Fig. 2

Change in evolutionary dynamics due to inclusion of intraspecies dynamics. When the fitness of Generous and Selfish types in both species is solely determined by the interactions that occur between species (in this case mutualism, $p=1$) then the dynamics are as reported previously [35]. Colours represent the initial states which result in an outcome favourable for species 1 (blue leading to (1S,2G)) and species 2 (red, leading to (1G,2S)). As in the same parameter configuration as the left panel ($p=1$) if the impact of intraspecies dynamics is reduced by a 1/3 ($p=0.666$), a qualitatively different picture is obtained. Two fixed points (solid black circles) are observed where both the Generous and Selfish types can co-exist in both the species. All initial state in the interior leads to either one of these fixed points (hence the lack of colours). However it is still possible to characterise the “successful" species. Depending on which of the two stable equilibria is reached, the species with fewer Generous individuals is the “successful" one. The null clines for species 1 (blue, horizontal) and for species 2 (red, vertical) are shown. The analysis was performed $d_1^{inter} = d_2^{inter} = d_1^{intra} = d_2^{intra} = 5$, $M_1 = M_2 = 1$, $b=2$, $c=1$ and $r_1 = r_2 /8$ for the interspecies mutualism game while additionally ${\tilde{b}}_1 = {\tilde{b}}_2 = 10$ and ${\tilde{c}}_1 = {\tilde{c}}_2 = 1$ and $\omega _1 = \omega _2 = 3/4$ for the two intraspecies games within each species. We chose a coexistence game to represent intraspecies interactions for both species denoted by $\rightarrow \bullet \leftarrow$. See Appendix C for different intraspecies interactions within each species and for different p values (Color figure online)

As a particular example, when the parameters used in Fig. 2 ($p=0.666$ panel) are held constant, there are two interior stable fixed points: both Generous and Selfish types can co-exist in populations, and at two different compositions. But external factors might mean that the coupling parameter p (the ratio of inter and intra species interactions) does not remain the same over time as in 6. In particular what happens when p changes in a continuous manner?

Fig. 3

Seasonal changes in the interspecies interactions affecting evolutionary dynamics within species. Interspecies interactions impact the fitness of the different types as in Eqs. 4, however instead of a static value for p, seasonality via a simple sine function as $p(t) = (1+\sin (at))/2$ is assumed. Here, a denotes how the seasonality time scale relates to the inter-intra-species interactions timescale. A large a denotes multiple bouts of mutualism affecting fitness for a given evolutionary time step while a small a denotes fewer such bouts within the same evolutionary time step. The trajectories shown in the panels are obtained by numerical integration with initial conditions $x_{1G} = x_{2G} = \{0.1,0.9\}$ and a step size of $\Delta x_{1G} = \Delta x_{2G} = 0.1$. The background colour is obtained by a finer grain of $\Delta x_{1\,G} = \Delta x_{2\,G} = 0.01$ and depicts the same outcomes as in Fig. 2, with gray representing the outcome where none of the edge equilibria are reached. For comparable or larger a the dynamics under oscillations can be captured by the average dynamics (at $p = 0.5$) however for small a a qualitatively different outcome is seen. Furthermore the phase in which the oscillating function begins is more important for smaller and smaller a especially if the stability of the fixed points changes as p changes (see Fig. 6 panel (b) x (b) across the p continuum)

Instead, consider a time-dependent function $p(t) = (1 + \sin (a t))/2$. The effect of such seasonality is demonstrated in Fig. 3. Introduction of seasonality maintains the two interior fixed points, but this is seen only when oscillations in which rates of change in p(t) are comparable, ($a=1$), or faster ($a=10$), with respect to the evolutionary timescale. For slower oscillations ($a=0.1$) cyclic behaviour emerges which is more prominent in species 2 than in species 1. Slow oscillations mean that the system spends longer close to the starting value of p and hence the initial phase of the p(t) oscillation becomes more important. This is especially interesting if the stability of the system is qualitatively affected. For example, if for a certain value of p the trajectories lose a certain type (reach an edge in the simplex), then a slow changing environment favouring the lost type would not be able to pull an evolutionary rescue.

2.5 Effect of Population Density on Interactions Dynamics

Population sizes change over time. Assuming that ecological changes are fast enough that they can be averaged, it is usually possible to ignore their effect on evolutionary dynamics. It is now possible to show that evolution can happen at fast timescales, comparable to those of the ecological dynamics [8, 37, 71, 79]. Hence it is necessary to consider not just evolutionary dynamics but eco-evolutionary dynamics together [17, 59].

Fig. 4

Population and evolutionary dynamics with combined inter-intra-species dynamics. As with the interactions described in Fig. 1 the two species consist of two types of individuals Generous and Selfish. Since the two species can in principle occupy different environmental niches, they can have non-overlapping population carrying capacities. The normalised carrying capacity in both species is 1 and $x_{1G} + x_{1S} + z_1 = 1$ (for species 1) where $x_{1G}$ and $x_{1S}$ are the densities of the Generous and Selfish types respectively. $z_1$ represents the remaining space into which the population of species 1 can still expand. For $z_1 = 0$ the species 1 is at its carrying capacity while for $z_1 = 0$ it is extinct. With the obvious substitutions the same equations apply for species 2

To include population dynamics in the previously considered scenario, $x_{1G}$ is reinterpreted as the fraction of Generous types and $x_{1S}$ as the fraction of Selfish types in species 1. Furthermore $z_1 = 1 - (x_{1\,G} + x_{1\,S})$ is the empty spaces in the niche occupied by species 1, and similarly for species 2 (Fig. 4). This approach has previously been explored in terms of social dilemmas in [33, 40], albeit for a single population. Here we adapt and modify it for two species. Hence the dynamics of the complete system is determined by the following set of differential equations,

$$\begin{aligned} {\dot{x}}_{1G}&= r_1 \; x_{1G} \; (z_1 \, f_{1G} \; - e_1) \end{aligned}$$

(7)

$$\begin{aligned} {\dot{x}}_{1S}&= r_1 \; x_{1S} \; (z_1 \, f_{1S} \; - e_1) \end{aligned}$$

(8)

$$\begin{aligned} {\dot{z}}_1&= - {\dot{x}}_{1G} \; - \, {\dot{x}}_{1S} \end{aligned}$$

(9)

for species 1, and similarly for species 2 with the proper index exchanges. Here $e_1$ and $e_2$ are the death rates of the two species. Setting $e_1 = z_1 {\bar{f}}_1$ and $e_2 = z_2 {\bar{f}}_2$ recovers the two species replicator dynamics as in Eqs. 6. In this scenario however fitnesses (and thus the average fitnesses) need to be re-evaluated as it becomes necessary to account for the presence of empty spaces. Thus now the fitnesses are as given in Appendix D. This is so because when forming a group of a certain size there may not be enough individuals to constitute it. Hence the effective group size can vary. Assuming that the minimal group size needs to be at-least 2 (to assure an interaction), it is allowed to vary up to the maximum number of players possible in a given game. Thus for example, the average fitness of the G strategy in species 1 is,

$$\begin{aligned} f_{1G}^{inter}({\textbf{x}}_{2})&= \sum _{s=2}^{d_1^{inter}-1} \left( {\begin{array}{c}d_1^{inter} -1\\ s-1\end{array}}\right) z_2 ^{d_1^{inter} -s} (1-z_2)^{s-1} P_G^{inter}(s,x_{2G},x_{2S},z_2) \end{aligned}$$

(10)

$$\begin{aligned} f_{1G}^{intra}({\textbf{x}}_{1})&= \sum _{s=2}^{d_1^{intra}-1} \left( {\begin{array}{c}d_1^{intra} -1\\ s-1\end{array}}\right) z_1 ^{d_1^{intra} -s} (1-z_1)^{s-1} P_G^{intra}(s,x_{1G},x_{1S},z_1) \end{aligned}$$

(11)

$$\begin{aligned} f_{1G}({\textbf{x}})&= p f_{1G}^{inter}({\textbf{x}}_{2}) + (1-p) f_{1G}^{intra}({\textbf{x}}_{1}) \end{aligned}$$

(12)

where the eventual group size s now depends on the empty space z. The payoff calculations now include an additional calculation, V(s, p, q, r) (see Appendix D). This function determines the probability of observing different group compositions whose size can be $s-1$, which can be less than the required game size due to the existence of empty spaces.

Since the population density is now variable $(x_{1G} + x_{1S})$ it is instructive to look at the proportion of Generous types within it over time. The dynamics can then be simplified by focusing on, $g_1 = x_{1G}/(x_{1G} + x_{1S})$ whose time evolution is given by,

$$\begin{aligned} {\dot{g}}_1&= r_1 z_1 g_1 (1-g_1) (f_{1G} - f_{1S}) \end{aligned}$$

(13)

$$\begin{aligned} {\dot{z}}_1&= r_1 (1-z_1) (e_1 - z_1 (g_1 f_{1G} + (1-g_1) f_{1S})) \end{aligned}$$

(14)

and similarly for ${\dot{g}}_2$ and ${\dot{z}}_2$ where everywhere we have $x_{1G} = g_1 (1-z_1)$ (with $x_{1S} = (1-g_1) (1-z_1)$) and $x_{2G} = g_2 (1-z_2)$ (with $x_{2S} = (1-g_2) (1-z_2)$) in the fitnesses as well.

Decomposing the model to focus on g and z, reflects the evolution (change in the proportion of generous individuals) and ecology (total population density/niche occupancy) of the system. Thus studying the dynamics of g and z together presents the eco-evolutionary dynamics of the system. Interactions at varying population densities affect group size formation which now includes the possibilities of player positions being left empty. Thus for smaller population densities the interaction groups are small and vice versa at higher densities. The effect of group size on evolutionary dynamics has been documented and can change the results qualitatively [66, 83]. Such a two species multi-type interaction system is a complicated, but likely, realistic depiction of most mutualisms observed in nature. Given this complexity, it is useful to consider the dynamics within the two species simultaneously.

To do this we take the most stable situation in which population dynamics are absent (Fig. 2) (defined by two internal stable equilibria) and incorporate population dynamics. The results are summarised in Fig. 5 where the evolutionary parameter (fraction of Generous in each species) is plotted against the ecological parameter (the number of the empty sites in the niche). The abundances of the two types, within the two species, change simultaneously. We therefore focus on the dynamics within a single species while the initial condition of the other species is kept fixed. In this manner it is possible to explore all initial conditions within the species of interest. For example, in the left panel of Fig. 5 the initial abundance of Generous types in species 2 is $g_2 = 0.5$ (starting at a density of $z_2 = 0.5$). In the phase space of ecology (density) and evolution (abundance), and starting at all possible initial conditions in species 1, coexistence is possible between the generous and selfish types depending on the impact of the intraspecies interactions (p) (Fig. 7). At the community level if two populations of species 1 end up in the two viable equilibria (as in Fig. 5, left panel) then there is a possibility for equilibrium selection. A similar process proceeds for species 2, as depicted in the right panel of Fig. 5 where all Generous can either disappear or take over the population. Thus while the intraspecies game would suggest coexistence between the Generous and Selfish types as defined within species, high impact of interspecies interactions can override the expectation (Fig. 7).

Fig. 5

The interplay of evolutionary strategies with population density for an intraspecies coexistence game with interspecies mutualism. With exactly the same parameters as that of Fig. 2, but with $p=0.9$ and symmetric death rates $e_1 = e_2 = 0.05$, numerically evaluated examples are presented. Equilibrium selection is possible for species 1 as two equilibria exist. Left Panel: shows the outcomes in species 1 when starting from 0.5 fraction of Generous individuals in species 2 at half carrying capacity $z_2 = 0.5$. While a proportion of the initial conditions lead to an extinction of the Generous types (turquoise) or very close to it (yellow), there exist another fixed point where Generous individuals take over the population at high density (blue). Right Panel: shows the dynamics in species 2 with 0.5 fraction of Generous individuals in species 1 with empty spaces in the proportion of $z_1 = 0.5$ as an initial condition. Generous individuals in species 2 also either go extinct or take over the population both at high densities. Indeed the evolutionary dynamics in the two species and the population dynamics as shown in the two panels needs to be ideally considered simultaneously. The numerical integration in both panels was run for 10000 steps and the endpoints assessed as equilibrium values. As the impact of intraspecies game (p) decreases, the equilibrium values do not allow for coexistence (see Fig. A.2) (Color figure online)

3 Discussion

The traditional view of mutualism is one in which there exists a harmonious relationship between species, each providing benefits to the other. However, classical theory on the evolution of cooperation demonstrates that such relationships are prone to exploitation. Such a situation can lead to ecological arms races, in which species play exploiter and exploited. While such reciprocal antagonistic interactions may persist in the short term, failure in the long term results in the loss of one or both species. Despite this fact, mutualisms appear to be remarkably common in nature. One possibility is that many putative mutualisms have been wrongly identified [28]; the other is that mutualisms are indeed common, but our understanding of the ecological and evolutionary factors shaping them is incomplete. A further possibility lies somewhere in between: mutualisms are common and persist over evolutionary time but are prone to failure over ecological time scales at any given locale.

Our investigation aligns with the latter possibilities and an expanding body of research that underscores the delicate equilibrium characterizing mutualistic relationships, where the balance between success and failure teeters precariously. It has become increasingly evident that secondary and even tertiary level interactions play a pivotal role in the long-term persistence of such associations [2, 37, 71, 79]. Our study and previous investigations highlight the limitations of models that depict mutualisms with linear species interactions. Understanding the mechanisms driving their maintenance becomes significantly challenging within such reductionist frameworks [19, 29]. However, models that incorporate even a modest degree of ecological and evolutionary realism, particularly nonlinear interactions within and between species, offer a more comprehensive perspective and demonstrate the potential for the persistence of mutualistic associations [2, 22]. As previously indicated in [50] and [77], it is reasonable to surmise that such intricate interactions extend to numerous other instances, further emphasizing their widespread prevalence and significance.

Incorporating intra- and inter-species interactions represents a significant stride in comprehending the intricate dynamics of eco-evolutionary processes. Nevertheless, it is essential to acknowledge that these factors alone might fall short of fully capturing the immense complexity inherent in such phenomena. The role of intraspecific interactions was already highlighted in [68]. Interspecific competition between two species of carrion flies could be alleviated by intraspecific aggression brought about by aggregation [45]. In contrast to competition, herein we focus of mutualistic interaction structures while making a similar argument about connecting intra and interspecific dynamics. We delve deeper into mutualistic associations and recognize the profound influence exerted by ecological fluctuations, capable of fundamentally reshaping the nature of species interactions [41]. To illustrate this point, we focus on an exemplary case study involving the hawk moth species Manduca sexta and the agave plant species Agave palmeri. Empirical investigations have illuminated the detailed intricacies of their mutualistic relationship, revealing the necessity of coordinated spatial and temporal behaviours. Intriguingly, these interactions exhibit density-dependent effects, hinting at underlying eco-evolutionary feedback mechanisms. Remarkably, the triadic involvement of a third species, Datura wrightii, further amplifies the complexity of this ecological interplay [73], a fact that we do not consider in our current model. Nevertheless, temporal coordination emerges as a critical factor, particularly in the context of pollination-related mutualisms [72]. Our findings, visually depicted in Fig. 3, underscore the paramount importance of seasonality in shaping the dynamics of these intricate interactions. However, it is imperative to recognize that the sensitivity of mutualistic associations to environmental conditions extends far beyond what might seem a mere exaggeration. As the weather systems of the globe shift due to anthropogenic activities, the phenological mismatch in pollination mutualisms are not the only ones reaching criticality. The ongoing climate change phenomenon has triggered a fundamental transformation in the coral-dinoflagellate symbiosis, steering it towards a precarious path of parasitism [7]. The repercussions of such mutualism breakdown are far-reaching and can potentially lead to catastrophic ecological meltdowns [18, 60, 92].

Similarly, the mutualism between Vibrio fischeri and the bobtailed squid appears to depend on a wide range of host and symbiont factors, although their specific contributions to the persistence of the mutualism are unclear. The squid provides a protected niche for V. fischeri, which in return provides light that the squid uses as a form of counter-illumination. A well studied ‘winnowing’ mechanism has been attributed [65, 82] to the separation of Vibrio from the rest of the microbiome and specially grown to high densities. Nevertheless, a closer analysis quickly leads to questions about how the interaction between the two species can be maintained in such a seemingly benign state. After all, V. fischeri colonises squid mucosal surfaces, which can be readily exploited even if the initial coloniser is a mutualist due to the rapid growth within the light organ. In addition, infections are typically established by a mixture of genotypes and the bacterium is horizontally transmitted [76]. These factors alone should drive V. fischeri to become virulent unless some checks are established [29, 52, 54]. Just how the squid avoids outright exploitation by the bacterium is not known, but the daily expulsion of bacteria by the squid results in diel changes in bacterial population density that may have more to do with limiting opportunity for within-host evolution than allowing the bacterial population to reach some optimal light-emitting status. What of population ecology? It is hard to imagine that V. fischeri mutants that take unfairly of host resources do not arise. Perhaps an opportunity for these types to persist is limited because of host density and limited opportunity for transmission. Thus inclusion of density-dependent dynamics are crucial when exploring the interactions between species Fig. 5.

As is evident throughout, our approach is implicitly multilevel. This allows community dynamics to be understood as a set of interactions between levels of selection [58, 87], and for within-species interactions to be understood in light of between-species interactions and vice versa. Such an approach shines light on a higher dimensional ecological space that is often overlooked: it facilitates predictions as to the range of parameters over which mutualisms can be maintained and allows exploration of the effects of community structure on the emergence of mutualisms in the absence of “game-changing" factors. This does not mean that game-changers, such as host sanctioning, are not important, but it shows that mutualisms can sometimes be understood simply in terms of community dynamics. From such an ecology-first perspective it becomes possible to understand how selective processes consequently shape the evolution of host sanctioning and various “lock and key" mechanisms that likely contribute to the long term persistence and refinement of mutualisms [88].

Finally, although we have included complex intra- and interspecies interactions along with density-dependence and phenology, the implementation outlined here has made several simplifying assumptions. For example, interaction terms used to define each player’s fitness are identical, and the threshold at which benefits are generated to each species are also identical, as are the number of players of each species. For simplicity, we have kept within-species strategies the same as that of between species. If there is no linkage between the traits of individuals when acting between or within species, then several interactions are possible. In a separate study, we have relaxed this assumption [91]. Assuming no linkage between strategies, we would have a two-population model with multiple strategies playing multiple games. Our previous work has shown that as long as the number of strategies is not more than two, the complex dynamics of multiple multiplayer games can be captured by a linear combination of the games. However, for multistrategy games, a dynamical inconsistency emerges [14, 15]. Thus the dynamics of the system can no longer be captured by a simple combination of the constituent interactions and demonstrates emergent phenomena.

We have explored all combinations of simple games in Appendix B and Fig. 6. With the theory developed, more complex scenarios can be readily analysed, including those involving addition of parameters to accommodate policing or sanctioning [16], and extension to finite populations [90]. Our framework is extensible to include a multi-species community context as discussed in the Introduction. That it is not necessary to invoke such complexity to explain the maintenance of mutualisms is an affirmation of the generality of our framework, but it also emphasises the role of the feedback between external ecological factors and evolutionary change. Since introduction of game theory in biology by [57], the primary focus of evolutionary game theory has been on antagonistic interactions and the resolution of conflicting interests to explain observations such as cooperation and mutualism. Numerous solutions have been proposed over the decades. Still, the explicit inclusion of ecological processes is emerging to be a necessity when aiming to provide a broader context [58, 86, 90]. Combining the dynamical nature of ecological processes with evolutionary games is thus the logical next step.

Appendices

Interspecies Evolutionary Dynamics

Traditional coevolutionary models consider interspecific dependence only [74, 75]. Since in our case each of the interactions between species are mutualistic and each species consists of two types of individuals Generous and Selfish, the following Snowdrift game is an appropriate representation of the interactions.

1.1 The Snowdrift Game

1.1.1 Two Player Setting

Two car drivers are stuck in a snowdrift. They must shovel away the snow (paying the cost c) to reach home (benefit b) but there are three possible outcomes to this scenario: One of the drivers shovels while the other stays warm in the car ($b-c$ and b); both the drivers share the workload and shovel away the snow ($b-c/2$ and $b-c/2$); neither of them gets out of the car and they both remain stuck (0 and 0).

Putting this game in perspective of the two species (i.e. the two drivers represent individuals of the two different species) we get the matrix,

(A.1)

where strategy G stands for being Generous and shoveling the snow while S stands for being Selfish and just sitting in the car. For $b=2$ and $c=1$ we recover the matrix used in [9].

For the snowdrift game in a single population for which the pairings are formed at random, there exists a single, stable internal equilibrium: the population will evolve to a polymorphism that is a combination of Generous and Selfish individuals. But in a two species system (pairs still random, but one from each species), this stable equilibrium turns into a saddle point: a small deviation from this fixed point leads the system to one of the stable fixed point where one of the species is completely Generous and the other one is completely Selfish [9] (as shown in Fig. 6, top).

Fig. 6

We assume that the interspecies game is set to a mutualistic interaction setting such as, $d_1^{inter} = d_2^{inter} = 5$, $b = 2$, $r_1 = r_2/8$, $M_1 = M_2 = 1$ and $c=1$. As for the intraspecies games we have $d_1^{intra} = d_2^{intra} = 5$ and ${\tilde{b}} = 10$ with (a) ${\tilde{c}} = 3$, $\omega = 3/4$, (b) ${\tilde{c}} = 1$, $\omega = 3/4$, (c) ${\tilde{c}} = 1$, $\omega = 4/3$ and (d) ${\tilde{c}} = 3$, $\omega = 4/3$, recovering the traditional scenarios of dominance, coexistence and bistability [39]

Fig. 7

Eco-evolutionary dynamics across two species. The interspecies game uses the same parameters as in Figs. 2 and 5 and the intraspecies dynamics allows for coexistence. Hence (b) scenario from Fig. 6. As the impact of interspecies interactions increases p, the proportion of initial conditions that result in coexistence equilibrium (yellow) (as evaluated after 10000 integrations steps) shrinks. Already after $p=0.9$ the Generous types either take over the population (blue), or go extinct (turquoise) (Color figure online)

1.2 Multiplayer Setting

Following Souza et al. [83], a multiplayer snowdrift game can be described by the payoff entries

$$\begin{aligned} \Pi _{1G} (k)&= {\left\{ \begin{array}{ll} b-\frac{c}{k} &{} \text {if } k \ge M \\ -\frac{c}{M} &{} \text {if } k < M \end{array}\right. } \end{aligned}$$

(A.2)

$$\begin{aligned} \Pi _{1S} (k)&= {\left\{ \begin{array}{ll} b &{} \text {if } k \ge M \\ 0 &{} \text {if } k < M. \end{array}\right. } \end{aligned}$$

(A.3)

All players gain a benefit b if the number of generous individuals in both species combined, k, is greater than or equal to the threshold M. For the generous individuals, their effort is subtracted from the payoffs. The effort is shared if the quorum size is met $\left( \frac{c}{M}\right)$, but is in vain for $k<M$. For two player games we had $M=1$ but multiplayer games provide the possibility of exploring this threshold aspect of collective action games. From these payoff entries we need to calculate the average fitnesses. For simplicity we just illustrate the fitnesses of the strategies in species 1. Species 1 plays a $d_1^{inter}$ player game. This means that every individual in species 1 interacts with $d_1^{inter} - 1$ individuals from species 2. Assuming random sampling the composition of the formed groups is given by a binomial distribution. Summing over all possible compositions of groups the average fitnesses of the two strategies in species 1 is obtained as,

$$\begin{aligned} f^{inter}_{1G} ({\textbf{x}}_2)&= \sum _{k=0}^{d_1^{inter} -1} \left( {\begin{array}{c}d_1^{inter} -1\\ k\end{array}}\right) x_{2G}^k (1-x_{2G})^{d_1^{inter} -1-k} \Pi _{1G}(k+1) \end{aligned}$$

(A.4)

$$\begin{aligned} f^{inter}_{1S} ({\textbf{x}}_2)&= \sum _{k=0}^{d_1^{inter} -1} \left( {\begin{array}{c}d_1^{inter} -1\\ k\end{array}}\right) x_{2G}^k (1-x_{2G})^{d_1^{inter} -1-k} \Pi _{1S}(k), \end{aligned}$$

(A.5)

and similarly $f_{2G}^{inter}$ and $f_{2S}^{inter}$ for species 2.

Note that here for the sake of notation we have assumed the same cost-benefit and threshold values for the two species. However along with the number of players $d_1^{inter}$ and $d_2^{inter}$, these parameters may be different for the two species. For asymmetric bi-matrix games there is a difference in the dynamics between the standard replicator dynamics and the alternative dynamics put forward by [56]. In the latter case the replicator equations cannot be simplified by removing the average fitness from the denominator and can give rise to qualitatively different dynamics. Then it is necessary to resort to difference rather than differential equations.

Intraspecies Evolutionary Dynamics

For elucidating intraspecies dynamics we focus on species 1 as the analysis is analogous for species 2. Within-species dynamics can in principle be completely different from the between-species interactions. The partitioning of individuals into two strategies follows the same partitioning as in inter-species interactions—Generous and Selfish. In principle it is possible to have two different labels for the strategies in the intraspecies interactions and the Generous and Selfish categories could be split among them. By this, the individuals are not constrained to follow the same strategy in both inter and intra species games. However for the sake of simplicity the same categorisation is assumed as in the inter-species level. A complete separation of the inter and intra strategy labels is implemented in another manuscript [91] where the ensuing complex dynamics are analysed in detail.

1.1 Synergy/Discounting Framework

Within species interactions are modelled using a general framework of costs and non-linear benefits [21, 39] which can encompass many different types of social interaction structures in a qualitative manner [63], i.e., dominance of either type, coexistence and bistability. Since the categorisation of the strategies at the intraspecies level is the same as that of the inter species level, for species 1 the frequencies of Generous and Selfish types are $x_{1G}$ and $x_{1S} = 1-x_{1G}$. The Generous and Selfish in species 1 play a $d_1^{intra}$ player game. Thus the fitnesses of the two types are defined as [39],

$$\begin{aligned} f^{intra}_{1G} ({\textbf{x}}_{1})&= \sum _{k=0}^{d_1^{intra} -1} \left( {\begin{array}{c}d_1^{intra} -1\\ k\end{array}}\right) x_{1G}^k (1-x_{1G})^{d_1^{intra} -1-k} \Gamma _{1G}(k+1) \end{aligned}$$

(A.6)

$$\begin{aligned} f^{intra}_{1S} ({\textbf{x}}_{1})&= \sum _{k=0}^{d_1^{intra} -1} \left( {\begin{array}{c}d_1^{intra} -1\\ k\end{array}}\right) x_{1G}^k (1-x_{1G})^{d_1^{intra} -1-k} \Gamma _{1S}(k). \end{aligned}$$

(A.7)

where the payoffs are given by,

$$\begin{aligned} \Gamma _{1S} (k) = \frac{{\tilde{b}}}{d_1^{intra}} \sum _{i=0}^{k-1} \omega ^i \end{aligned}$$

(A.8)

$$\begin{aligned} \Gamma _{1G} (k) = \Gamma _{1S} (k) - {\tilde{c}}. \end{aligned}$$

(A.9)

Thus Selfish types gain a fraction of the benefit which is scaled by the factor $\omega$, which determines whether the benefits are linearly accumulating ($\omega =1$) for increasing number of Generous individuals, synergistically enhanced ($\omega >1$), or saturating ($\omega <1$). Note that the costs and benefits of the within-species game need not be the same as those between species ($b\ne {\tilde{b}}$ and $c \ne {\tilde{c}}$).

Combined Evolutionary Dynamics

The average payoffs are assumed to be a linear combination of the interspecies and intraspecies interactions where the parameter p determines the strength of each of the interactions such that,

$$\begin{aligned} f_{1G} ({\textbf{x}})&= p f^{inter}_{1G} ({\textbf{x}}_2) \; + \; (1-p) f^{intra}_{1G} ({\textbf{x}}_1) \end{aligned}$$

(A.10)

$$\begin{aligned} f_{1S} ({\textbf{x}})&= p f^{inter}_{1S} ({\textbf{x}}_2) \; + \; (1-p) f^{intra}_{1S} ({\textbf{x}}_1). \end{aligned}$$

(A.11)

Following the same procedure for the two strategies in species 2 leads to the average fitness

$$\begin{aligned} {\bar{f}}_1 ({\textbf{x}})&= x_{1G} \; f_{1G} \; ({\textbf{x}}) + x_{1S} f_{1S}({\textbf{x}}) \end{aligned}$$

(A.12)

$$\begin{aligned} {\bar{f}}_2 ({\textbf{x}})&= x_{2G} \; f_{2G} \; ({\textbf{x}}) + x_{2S} f_{2S}({\textbf{x}}). \end{aligned}$$

(A.13)

The time evolution of Generous types in both the species provides complete dynamics of the system. However since the two interaction species are by definition different organisms, they can have different rates of evolution. Thus if species 1 evolves at the rate $r_1$ while species 2 at rate $r_2$ then we have,

$$\begin{aligned} {\dot{x}}_{1G}&= r_1 \; x_{1G} \; \left( f_{1G}({\textbf{x}}) - {\bar{f}}_1({\textbf{x}}) \right) \end{aligned}$$

(A.14)

$$\begin{aligned} {\dot{x}}_{2G}&= r_2 \; x_{2G} \; \left( f_{2G}({\textbf{x}}) - {\bar{f}}_2({\textbf{x}}) \right) . \end{aligned}$$

(A.15)

Population Dynamics

For brevity we begin with the description of population dynamics in species 1. The two types in species 1, Generous and Selfish need not sum to 1, i.e., the population may not always be at carrying capacity. Hence if the empty space in the niche occupied by species 1 is $z_1$, then we have $x_{1G} + x_{1S} + z_1 = 1$ where $x_{1G}$ and $x_{1S}$ are the densities of Generous and Selfish types. In species 1 the population dynamics is determined by,

$$\begin{aligned} {\dot{x}}_{1G}&= r_1\; x_{1G} (z_1 f_{1G}({\textbf{x}}) - e_1({\textbf{x}})) \end{aligned}$$

(A.16)

$$\begin{aligned} {\dot{x}}_{1S}&= r_1\; x_{1S} (z_1 f_{1S}({\textbf{x}}) - e_1({\textbf{x}})) \end{aligned}$$

(A.17)

$$\begin{aligned} \dot{z_1}&= - {\dot{x}}_{1G} - {\dot{x}}_{1S} \end{aligned}$$

(A.18)

and similarly in species 2, where $e_1$ and $e_2$ are the death rates for the two species. For the special case of $e_1({\textbf{x}}) = z_1 \bar{f_1}({\textbf{x}})$ where $\bar{f_1}({\textbf{x}}) = \frac{x_{1\,G} f_{1\,G} + x_{1\,S} f_{1\,S}}{x_{1\,G} + x_{1\,S}}$ is the average fitness of Species 1 (and similarly for $e_2$) the two species reduce to replicator dynamics as in Eqs. A.15. The fitnesses however need to be re-evaluated in this scenario. For example in species 1 the fitness for type G is,

$$\begin{aligned} f_{1G}^{inter}({\textbf{x}}_{2})&= \sum _{s=2}^{d_1^{inter}-1} \left( {\begin{array}{c}d_1^{inter} -1\\ s-1\end{array}}\right) z_2 ^{d_1^{inter} -s} (1-z_2)^{s-1} P_G^{inter}(s,x_{2G},x_{2S},z_2) \end{aligned}$$

(A.19)

$$\begin{aligned} f_{1G}^{intra}({\textbf{x}}_{1})&= \sum _{s=2}^{d_1^{intra}-1} \left( {\begin{array}{c}d_1^{intra} -1\\ s-1\end{array}}\right) z_1 ^{d_1^{intra} -s} (1-z_1)^{s-1} P_G^{intra}(s,x_{1G},x_{1S},z_1) \end{aligned}$$

(A.20)

$$\begin{aligned} f_{1G}({\textbf{x}})&= p f_{1G}^{inter}({\textbf{x}}_{2}) + (1-p) f_{1G}^{intra}({\textbf{x}}_{1}) \end{aligned}$$

(A.21)

and similarly for type 1S where the payoff functions are defined as,

$$\begin{aligned} P_G^{inter}(s,p,q,r)&= \sum _{k=0}^{s-1} V(s,p,q,r) \Pi _{1G}(k+1) \end{aligned}$$

(A.22)

$$\begin{aligned} P_G^{intra}(s,p,q,r)&= \sum _{k=0}^{s-1} V(s,p,q,r) \Gamma _{1G}(k+1) \end{aligned}$$

(A.23)

$$\begin{aligned} P_S^{inter}(s,p,q,r)&= \sum _{k=0}^{s-1} V(s,p,q,r) \Pi _{1S}(k) \end{aligned}$$

(A.24)

$$\begin{aligned} P_S^{intra}(s,p,q,r)&= \sum _{k=0}^{s-1} V(s,p,q,r) \Gamma _{1S}(k) \end{aligned}$$

(A.25)

where $V(s,p,q,r) = \left( {\begin{array}{c}s-1\\ k\end{array}}\right) \left( \frac{p}{1-r}\right) ^k \left( \frac{q}{1-r}\right) ^{s-1-k}$ is the probability of having a k Generous individuals and $s-1-k$ Selfish individuals in the inter(intra) species game. and the actual payoffs are calculated as per Eqs (A.3) and (A.9).