Appendix A: Individual selection in ecosystems
A mutation to an individual in species i that decreases the competitive effect, ω
ij
, of species j on species i directly affects the fitness of the individual carrying the mutation and not other individuals in species i, and can thus be favoured by individual selection. It is only changes to traits that directly affect the growth rate of an individual compared to the rest of the individuals in the species can be affected by individual-level selection. Traits that increase the growth rate of all individuals in the species equally have no differential individual benefit (despite conferring benefit to the species as a whole). In particular, a mutation to a trait in an individual in species i that changes its competitive effect, ω
ji
, on some other species j, e.g. decreasing the density of a competitor species, may thereby indirectly increase the growth rate of species i. But this will benefit all individuals in species i, not just the mutant, and therefore has no differential selective benefit to the individual that bears the mutation [77]. Likewise, the competitive effect of species j on species i may, by virtue of normalising ecological constraints, be decreased as a side-effect of increasing the competitive effect of species k on species i. But again this would not be favoured by individual selection as the benefit is felt by all individuals in species i (conversely, changes to ω
ij
could be selected under individual selection even though, as a result of indirect effects through changes in density of other species or through normalising ecological constraints, their net effect is to decrease the density of their own species). It is therefore only direct effects on individual fitness that are taken into account by the selection coefficient described here; i.e., Eq. 2. evaluates the change in growth rate of individuals in species i due to changes in ω
ij
and not ω
ji
, and furthermore, only changes to ω
ij
caused by positive selection coefficients, not those caused by indirect normalisation effects. This correctly disregards any changes to a species growth rate that occurs as an indirect side-effect of altering the density of some other species.
Appendix B: The relationship between rate of adaptation and product of species densities in more complex cases
In the main text, the rate of adaptation, v
ij
, of each interspecific interaction coefficient is modelled with Eq. 3 corresponding to the case where there is no interference between simultaneously segregating alleles at different loci. In large sexual populations with linked loci, the rate of adaptation will depend on the type of recombination, recombination rate, population size, the mutation rate and magnitude of mutations. Here we compare the rate of adaptation of an interaction coefficient for three different models. In each case, the rate of adaptation, v
ij
, of an interspecific interaction coefficient describing the fitness effect of species j on species i, is \(v_{\textit {ij}} = x_{i} \mu \bar {P}\), where x
i
μ is the rate with which beneficial mutations arise in species i, and \(\bar {P}\) is the average probability that a single new mutation will ultimately fix (see main text). In all cases, \(\bar {P}\) is a function of the selection coefficient \(s_{i} = \frac {m_{i}}{ k_{\textit {ie}}} g x_{j}\) (Eq. 2, main text) where m
i
is the intrinsic net growth rate of species i, k
i
e is the carrying capacity of species i in environment e, and g is the change in the interaction coefficient due to an individual mutation. Here we write s
i
=β
x
j
, for clarity of the comparisons that follow.
Case a) No interference
In simple cases when there is no interference between simultaneously segregating alleles at different loci (e.g. where genes are under weak selection per locus, free recombination and the linkage disequilibria among alleles sweeping to fixation are negligible), the probability of fixation, \(\bar {P} =s_{i}\). Thus, as per Eq. 3 main text:
$$ v_{ij} = \beta \mu x_{i} x_{j} $$
((4))
where x
i
is the density of species i, x
j
is the density of species j and μ is the beneficial mutation rate.
Case b) Linked genes on a linear genome
Weissman & Barton [79] consider the effects of interference between linked genes on a linear genome. Here the genomic rate of fixation of beneficial mutations is ([79] Eq. 7):
$$ v = \frac{v_{0}}{1+2 v_{0}/R} $$
((5))
where, v
0 is the genomic rate of fixation of beneficial mutations in the absence of interference and R is the total genetic map length in Morgans. The authors use the approximation v
0=2x
μ
s, where x is species density and s is the selection coefficient. With s
i
=β
x
j
as before, this gives the rate of adaptation on an interaction:
$$ v_{ij} = \frac{2 \beta \mu x_{i} x_{j}}{1 + 4 \beta \mu x_{i} x_{j} /R} $$
((6))
Case c) Occasional outcrossing
Neher et al. [78] study the rate of adaptation in unlinked loci in facultative sexuals where the rate of outcrossing is very small. Whereas Weissman and Barton examine the case of obligately sexual populations, this case represents occasionally/facultatively sexual populations (e.g. plants). On condition that r
2/s
2≫4x
μ, the rate of accumulation of beneficial mutations in this case is given by ([78] Eq. 12b):
$$ v \approx x \mu s^{2} \left(1 - \frac{4 x \mu s^{2}}{r^{2}} \right) $$
((7))
where r is the outcrossing rate. With s
i
=β
x
j
as before, this gives the rate of adaptation on an interaction:
$$ v_{ij} \approx x_{i} \mu (\beta x_{j})^{2} \left(1 - \frac{4 x_{i} \mu (\beta x_{j})^{2}}{r^{2}} \right) $$
((8))
Comparison of the three cases
Figure 7 plots the rate of adaptation v
ij
as a function of x
i
and x
j
for these three different cases. We observe that case a, where the rate of adaptation is directly proportional to the product x
i
x
j
as modelled in our simulations, and the two more complex cases (b and c) are all qualitatively similar. Although in some cases the absolute rate of adaptation is more strongly limited by the recombination rate than the mutation supply or the strength of selection, for example [79], the relative rates of adaptation are still determined largely by the product of x
i
and x
j
. More specifically, all three cases have the essential characteristic that the rate of adaptation is zero when either x
i
or x
j
is zero, and otherwise, the greater the value of one, the greater the rate of increase with the other. Thus, although the shape of the alternate functions differs from ours, the essential behaviour is preserved. Intuitively, mutations must be both created and selected for an interaction coefficient to evolve.
Appendix C: Additional methods: Normalisation, variable environments, measuring ecological attractors and model parameters
Normalisation methods
In each evolutionary step, all interaction terms in Ω(t) are updated by natural selection according to Eq. 3 to produce Ω
′(t) and then renormalised to produce Ω(t+1). Renormalisation preserves the conditions that for each species i and all other species \(j (j \ne i),\sum _{j=1}^{N} \omega _{\textit {ij}} (t) = Q_{i}\), and \(\sum _{j=1}^{N} \omega _{\textit {ji}} (t) = Q_{i}\), where Q
i
<0 is a constant for each species. Specifically, an iterative row and column normalisation (below) is applied to M(k=0)=Ω
′, until the values of M converge within a specified accuracy, i.e. \((\sum _{\textit {ij}} (m_{\textit {ij}} (k+1) - m_{\textit {ij}} (k))^{2} < 10^{-5}\), where k is the iteration counter, as follows:
$$ M(k+1)=\text{column\_norm} (\text{row\_norm} (M(k))) $$
((9))
where row_norm(m
ii
) = m
ii
, column_norm(m
ii
) = m
ii
, i.e. self-interactions are unaffected, and
$$ \text{row\_norm}(m_{ij(i \ne j}) = \frac{m_{ij}(k)}{\sum_{j=1 (j \ne i)}^{N} m_{ij} (k)} $$
((10))
and
$$ \text{column\_norm}(m_{ij(i \ne j}) = \frac{m_{ij}(k)}{\sum_{i=1 (i \ne j)}^{N} m_{ij} (k)} $$
((11))
Variable environments
We investigate the effect of variable environments as follows. The carrying capacity of the i
th species in a default ecological environment, E
0, is k
i0. For simplicity in our simulations we let k
i0=k
0, for all i, where k
0 is a constant. Prior to the evolution of interactions, this causes all species to equilibrate at the same density. To model the evolution of an ecosystem under varying environmental conditions that force or drive the ecosystem to adopt different ecological states, we define two other environmental conditions that alter carrying capacities. The pattern of equilibrium species densities under one environmental condition, E
1, increases the carrying capacity of some species to k
0+α and decreases others to k
0−α, where α=0.1. In E
2, a different subset of species is increased/decreased in a similar manner. See Fig. 2 main text.
Measuring ecological attractors
We examine the ecological attractors in the ecosystem by Monte Carlo sampling, i.e., by repeatedly setting the species densities to random initial conditions and running to an equilibrium. To measure the inherent attractors induced by evolutionary changes, this sampling is carried out in the absence of environmental forcing – i.e., in E
0. In some experiments we also investigate the amount of environmental forcing required to push the ecosystem out of equilibrium in one pattern of species densities and into the attractor basin of another stable equilibrium. Whenever, as here, interactions control the correlation of species densities and not their absolute densities, the complement of any attractor pattern is also necessarily an attractor [57, 58, 67]. However, so long as initial conditions are more similar to the past states experienced during evolution than the opposite of those past states these unnatural attractors are precluded. Accordingly, we examine initial conditions, x, satisfying the condition (|x−E
1|<|x−E1′|) and (|x−E
2|<|x−E2′|) where E
′ is the inverse of E (i.e. E
′ = \(2 \bar {E}-E)\).
Model parameters
N=400, number of species. m
i
=0.5, growth rate of all species. s(t=0)=0.1, initial species densities. k
0=10, a parameter governing the extrinsic component of carrying capacity in E
0. α=0.1 increment/decrement of particular carrying capacities in environments E
1 and E
2. T=1, number of evolutionary changes applied in each environment before switching. τ=5000, number of ecological timesteps (Eq. 1) between ‘initial’ and ‘final’. g=0.01, constant of proportionality in selection-limited evolution (Eq. 3)
Interaction coefficients are initialised as follows:
$$\omega_{ij}(t=0) =\left\{ \begin{array}{ll} -1, & \text{if}\, i=j (\text{i.e. self interactions})\\ -0.2, & \text{otherwise} \end{array} \right. $$
\(Q_{i} = \sum _{j=1 (j \ne i)}^{N} \omega _{\textit {ij}}(t=0)\), normalisation constant (the sum of the non-self interactions in any one row/column remains equal to their sum at time t=0).
The quantitative values of these parameters will naturally have quantitative effects on the behaviour of the eco-evolutionary dynamics that we simulate. Since the simulations are a phenomenological model of ecosystem evolution, what matters is the relative rather than absolute rates of adaptation on different interaction coefficients – in particular, which interactions increase, which decrease and which remain largely unchanged. This pattern, and its sensitivity to different modelling choices, is investigated in Appendix D.
Appendix D: Equivalence of Hebbian and evolved changes in more complex cases
In the main text the rate of adaptation of each interspecific interaction coefficient is modelled with Eq. 3 corresponding to the case where there is no interference between simultaneously segregating alleles at different loci. Appendix B shows that the characteristics of the rate of adaptation in more complex cases is qualitatively similar although they are quantitatively different. Here we simulate evolution using these alternative models and incorporating normalising ecological constraints. Figure 8 shows that the quantitative differences in the three equations do not alter the pattern of positive, negative and neutral changes that are produced in the evolving interaction matrix. Specifically, the pattern of changes in interactions have the same direction as the Hebbian model in all cases. Accordingly, there will be parameter ranges where they produce the same distributed memory phenomena in the ecosystem. Investigations of quantitative differences remain for future work.
Appendix E: Response to environmental forcing that is not similar to environments experienced during evolution
Figure 9 shows that an ecosystem can exhibit a non-catastrophic response when forced in arbitrary directions (b) and simultaneously exhibit hysteresis and catastrophic regime shifts when forced in directions that have been experienced previously over evolutionary time (a). This emphasises that the evolved ecological memory causing the switching behaviour is conditioned by the systems’ evolutionary history, and thus causes recall (or recognition) of a specific point in a multi-dimensional space of species densities, rather than a general stability/instability property resulting from unorganised or arbitrary evolutionary changes.
Appendix F: Development and breakdown of multiple attractors over long evolutionary timescales
Figure 10 shows how the attractors of the ecosystem change over evolutionary time in Experiment 2. Interestingly, we see that in the long term the two-attractor state is unstable because, rather than reinforcing the ecological patterns that are ‘forced’ by the external environment, the system begins to reinforce its own patterns of behaviour [58], and positive feedback causes one (slightly stronger) attractor to outcompete the other (Fig. 11).
Appendix G: Empirical tests for distributed learning in ecosystems
The dynamical behaviours we observe in the evolved ecosystem are consistent with ecological memory, alternate ecological states, succession dynamics, assembly rules, regime changes and founder effects observed in natural ecosystems. These behaviours follow from simple component principles (i.e. the availability of heritable variation in inter-specific interactions, and the presence of ecological constraints or evolutionary trade-offs) and direct evidence for these behaviours is testable. For example, consider the evolution of a small microbial community. Given a culturable community with stable coexistence dynamics, we could first test whether it has i) one or ii) alternative stable states. This requires sampling many different initial species compositions and allowing species densities to equilibrate. i) If a single state, we can then force the system into a different state (‘alternate ecosystem state’, [17]) – e.g. by changing temperature, nutrient influx – and hold it there for evolutionary time. Then remove the forcing and retest for multiple attractors (‘alternate community states’). If a memory has been conditioned by this forcing then a new attractor will be exhibited. ii) If the system initially has more than one attractor state, then we can estimate the basin size for each attractor by counting the number of different initial conditions that arrive at one or the other. By leaving the system in one attractor over evolutionary time this should increase the relative basin size in proportion to the time spent in that attractor. Next we need to assess the extent to which such a memory is collective or merely the sum of individual memories. This can be done by swapping-in evolved species for species in the original community one-by-one and assessing the relative contribution of individual and collective genetic changes on the dynamical behaviour of the system.
Appendix H: Asymmetric interactions, the importance of normalising ecological constraints, and other future work
One important aspect of evo-eco dynamics that is highlighted by this model is the importance of normalising ecological constraints or evolutionary trade-offs for collective behaviours. These constraints prevent a species A from benefiting from the presence of species B without also becoming dependent on B. That is, it is not just the case that A grows faster in the presence of B, but that A’s growth is slower when B is absent. Under these conditions, changes to interactions do not merely increase the growth of each species in a manner that is sensitive to its ecological context, but more specifically, they modify correlations between species densities. We assume in the present model that an adaptation that, for example, decreases the niche overlap with one species increasing the niche overlap with others. But the extent to which species evolve dependencies rather than just (context-sensitive) individual advantages in natural ecosystems is an empirical matter – and from this work we recognise it as a matter that is centrally important to the possibility of collective behaviours that are more than the sum of the individual behaviours.
This paper has investigated only competitive interactions and has not investigated mutualistic interactions or asymmetric interactions such as characteristic of trophic, e.g. predator-prey, relationships. The observation that selected changes to interactions are Hebbian does not depend on them being symmetric (or competitive). That is, Eq. 3 is not sensitive to any assumptions about the initial values of interaction coefficients, e.g., whether ω
ij
and ω
ji
are equal or even have the same sign, and therefore applies to predator-prey relationships as well as symmetric competitive interactions. Eq. 3 also shows that the selective pressures on changes to interactions are symmetric i.e., Δ
ω
ij
=Δ
ω
ji
(except for the influence of individually-varying carrying capacities), so there is no systematic reason for interactions to become asymmetric over evolutionary time. In the examples investigated in this paper the interaction coefficients are initialised symmetrically and, accordingly, they remain approximately symmetric. The evolutionary model could be applied to asymmetric interactions, but asymmetric interactions introduce the possibility of non-fixed point attractors, e.g. cycles, that complicate the behaviour of the eco-evolutionary dynamics and their measurement considerably. (We note that where ω
ij
and ω
ji
differ, the addition of multiple symmetric changes through natural selection will make them less asymmetric over evolutionary time, i.e., bring the ratio of these terms closer to 1, and could evolve them to take the same sign even when they started out with opposite signs. This implies that the effect of evolutionary change would be to increase the stability of the ecological dynamics and reduce or remove chaotic or cyclic attractors over time).
We have assumed that each interaction coefficient is independently modifiable whereas in natural populations traits may affect many interactions simultaneously. Here we chose to investigate scenarios where none of the interaction coefficients reach zero or go positive (which is possible in principle despite the normalisation employed). The equations used exhibit unstable behaviour in this case and a different approach to modelling would be required to handle mutualistic interactions. In natural populations one member of a population can gain selective advantage by changing its relationship to other members of its own species, but our simulations have fixed self-interactions at –1 and have investigated only the evolution of interactions with members of other species.
A key technical distinction between the recent work on associative memory in gene networks [67] and the models utilised here is that the Lotka-Volterra equations represent unsigned (positive) state variables, as is natural for species densities, rather than signed (positive and negative) state variables representing under- or over- expressed gene activity (compared to some normal level). Although it is possible and common to model interesting dynamical behaviours using either signed or unsigned state variables in neural networks, the use of unsigned variables means that Hebb’s rule, or natural selection, will only alter interactions in one direction, i.e., the product x
i
x
j
is always positive (although crucially it may have different magnitudes). The assumption of normalising constraints that cause some interactions to become more competitive as a side effect of others becoming less competitive is thus important to the results that we have shown.
In particular, as mentioned above, the assumption of these normalising constraints means that changes to interactions, although motivated by increases in individual growth rates, have the effect of (also) altering the dependency of one species on another. Without these constraints, the effect of unconstrained changes to interactions is to make high density species fitter in all conditions, rather than making them dependent on the simultaneous high density state of specific species (and hence less fit in some conditions). It is therefore important for future work to investigate how different ways of modelling such constraints impact the behaviours illustrated here. For example, rather than a Lotka-Volterra model, a stoichiometric model of species interactions may alleviate the need for an explicit normalisation mechanism.
Assuming that ecological dynamics (i.e., changes in species density) are much more rapid than evolutionary changes (i.e., genetic changes affecting the coefficients of inter-species fitness dependencies) [91], most evolution occurs whilst ecological dynamics are at or near equilibrium, as modelled here. The behaviour of evo-eco dynamics when these processes have more similar timescales [35] deserves attention. However, the fact that we model varying ecological conditions, causing the ecosystem to visit more than one ecological equilibrium, means that the interaction of ecological and evolutionary dynamics is non-trivial even though their timescales are kept separate in our simulations (following [24]). Moreover, any model assuming a single ecological attractor will overlook the interesting behaviours modelled here, regardless of whether the timescales are separated or similar.