# Prevalence of mutualism in a simple model of microbial co-evolution

Luciano Stucchi  Universidad del Pacífico, Lima, Peru  Group of Complex
Systems, Universidad Politécnica de Madrid, Madrid, Spain  Javier Galeano
Juan Manuel Pastor  Group of Complex Systems, Universidad Politécnica de
Madrid, Madrid, Spain  Jose María Iriondo  Biodiversity and Conservation Area,
ESCET, Universidad Rey Juan Carlos, Madrid, Spain  José A. Cuesta  Grupo
Interdisciplinar de Sistemas Complejos (GISC), Madrid, Spain  Dept.
Mathematics, Universidad Carlos III de Madrid, Leganés, Madrid, Spain
Instituto de Biocomputación y Física de Sistemas Complejos (BIFI),
Universidad de Zaragoza, Zaragoza, Spain

###### Abstract

Evolutionary transitions among ecological interactions are widely known,
although their detailed dynamics remain absent for most population models.
Adaptive dynamics has been used to illustrate how the parameters of population
models might shift through evolution, but within an ecological regime. Here we
use adaptive dynamics combined with a generalised logistic model of population
dynamics to show that transitions of ecological interactions might appear as a
consequence of evolution. To this purpose we introduce a two-microbial toy
model in which population parameters are determined by a bookkeeping of
resources taken from (and excreted to) the environment, as well as from the
byproducts of the other species. Despite its simplicity, this model exhibits
all kinds of potential ecological transitions, some of which resemble those
found in nature. Overall, the model shows a clear trend toward the emergence
of mutualism.

††preprint: APS/123-QED

##  I Introduction

In 1966, Jeon and Lorch were conducting experiments with a population of
Amoeba proteus. One of the strains got infected by X-bacteria, a Gram-
negative, rod-shaped bacteria related to _Legionella_ sp. and _Pseudomonas_
sp. (for which the name _Candidatus Legionella jeonii_ sp. nov. was later
proposed [1]). The few survivors to the infection retained the bacteria as
parasites [2]. Over the course of a few generations though, they became
endosymbionts, providing amoebas heatshock protection [3, 4].

Arguably, this is one of the most spectacular cases of an evolutionary
transition in an ecosystem from an antagonistic interaction to an obligate
mutualism, mainly because it was witnessed in real time. However, other
similar transitions are well documented in the literature from different
sources. For instance, phylogenetic data collected from one of the best-known
mutualistic systems—plants and pollinators—dating back 90 million years,
allowed [5] to show that the wasp Ceratosolen galili, which belongs to a genus
of active pollinating species, not only does not pollinate anymore but has
become a parasitic species of figs.

Likewise, performing phylogenetic analysis on 15 species of aphids of the
genus Chaitophorus, [6] concluded that the relationship between ants and
aphids may have changed at least five times during the history of their life,
thus going through all shades between mutualism and antagonism. There is no
consensus whatsoever on whether the ancestral mutualistic relationship between
these species was facultative or obligatory, but it seems reasonable to think
that, before these species developed any mutualistic relationship, ants
predated on aphids because this behaviour is still observed in all ant species
in the appropriate environmental conditions [7, 8, 9].

The above-mentioned examples are only a few well-studied ecological systems
that show transitions between different types of ecological interactions. In
all of them, the evolutionary nature of these transitions has been established
either through experiments or using phylogenetic analyses [10], but their
presence suggests that many more may have occurred in nature. This fact raises
many questions: Are these evolutionary transitions a common phenomenon or a
rarity? Why do ecological interactions move in one direction or another? Would
it be possible to predict in which direction a type of interaction will move,
if at all?

One way to explore the answers to some of these questions is to design
mathematical models that include two different time scales: a short one that
accounts for the usual population dynamics, and a long one that includes
Darwinian evolution. The latter can be accomplished by applying adaptive
dynamics to the parameters of the population model [11, 12, 13]. However, the
former is more problematic because standard population models deal with
different types of interactions in very different ways.

When it comes to modelling competition and predation, Lotka-Volterra
population equations are a very convenient choice [14, 15, 16]. Despite their
simplicity, these equations yield a rich set of biological predictions, to the
point that they can be used as a prototype for more realistic models as well
as a tool for interpreting complex observations. However, the Lotka-Volterra
model meets serious difficulties to accommodate mutualism because positive
interactions between species induce spurious feed-back loops that may drive
populations out of control [17]. The addition of Holling-type functional
responses [18] is the usual way to control an unbounded growth of the
populations, at the expense of rendering the model analytically intractable.
Recently though, a new population model has been introduced in which both,
species’ intrinsic growth and interspecific interactions, are limited via
logistic terms [19, 20]. The resulting equations are amenable to analytic
treatment and, at the same time, provide sensible results whether the
interactions are of antagonistic or mutualistic nature. These two nice
features make this model particularly suitable to study the evolution of
ecological interactions.

Variation in ecological interactions may take place as a result of multiple
factors [21]. For instance, the outcomes of ecological interactions may depend
on the age and/or life cycle stage of the individuals of the interacting
species. Similarly, phenotypic differences between the individuals of each of
the interacting species, may alter the impact of one species on the other.
These differences can result from genotypic diversity in the interacting
populations and/or phenotypic plasticity. Environmental conditions might also
determine how members of two populations interact with each other, allowing
transitions from mutualism to competition or even extinction [22]. In order to
focus on the general patterns derived from population dynamics and
evolutionary adaptation, we have chosen a basic microbial system in which
there are no differences among the individuals of each of the interacting
species as a result of age or life cycle stage. Furthermore, there is one
single genotype per interacting species and the environmental conditions are
fixed. Therefore, all individuals of each interacting species have the same
phenotype. Also for simplicity we will assume that mutations are so rare that
they either get quickly fixed in the whole population or disappear before a
new one occurs. In this regime, a suitable theoretical tool to study evolution
is adaptive dynamics [11, 12, 13].

Accordingly, our toy model consists of two microbial species that interact via
the resources they consume and excrete. They compete for resources from the
environment—when both of them use the same resource—but they can also cross-
feed from resources excreted by the other species. This interplay between
resource consumption and excretion will provide specific functional forms for
the parameters of Stucchi  _et al._ ’s population equations, automatically
yielding natural trade-offs between these parameters [23]. Adaptive dynamics
will then introduce the Darwinian mechanism for the evolution of the
parameters of the model.

Simplistic as it may be, this model provides sufficient complexity and
flexibility as to show a great diversity of evolutionary transitions. It
reveals a global trend toward the appearance of mutualistic interactions from
initially competitive or antagonistic scenarios, and exhibits evolutionary
pathways akin to some documented in the literature.

In what follows we provide a detailed account of the model for two species,
starting with a brief description of Stucchi  _et al._ ’s ecological model and
further connecting its phenomenological parameters with the microscopic
interactions between the species. The equations of the adaptive dynamic for
this model are developed in the Appendix B. We end by providing and analysing
the results of extensive numerical simulations.

##  II A toy model for two interacting microbial species

We will consider two microbial species whose populations evolve in time
according to the _generalised logistic model_ [20], briefly described in Eq.
(13). This model is suitable for our purpose because it allows for all kinds
of ecological interactions (whether beneficial or detrimental) and, at the
same time, populations are always limited by carrying capacities. Generalized
Lotka-Volterra equations lack this property, so that they end up having
difficulties in handling mutualism—hence introducing a spurious bias toward
competing interactions.

If we specialize this model for just two species, and ignore intraspecific
cooperation or direct competition (i.e., $b_{ii}=0$), the two equations that
describe this minimal ecological community are

$\begin{split}\dot{N_{1}}&=N_{1}(r_{1}-a_{1}N_{1})+N_{1}(1-c_{1}N_{1})b_{12}N_{2},\\\
\dot{N_{2}}&=N_{2}(r_{2}-a_{2}N_{2})+N_{2}(1-c_{2}N_{2})b_{21}N_{1}.\end{split}$
(1)

where the constants $r_{i}$ are the intrinsic growth rates, the $a_{i}$
account for the effect of intraspecific competitions, the $b_{ij}$ are the
interspecific interaction coefficients, and the $c_{i}$ are coefficients to
saturate the effect of interspecific interactions. This system of differential
equations is able to describe every kind of ecological interaction. Table 1
shows the kind of relationships of the species with the environment (through
the signs of the intrinsic growth rates $r_{i}$) and between them (through the
signs of the interaction constants $b_{ij}$) which correspond to the most
usual ecological scenarios. Depending on the signs of $b_{ij}$ the interaction
between both species can be _mutualistic_ (both positive), _antagonistic_ (one
positive and one negative, also called _predation_ or _parasitism,_ depending
on the context), or _competitive_ (both negative). On the other hand, the
relation with the environment determines whether mutualism is _facultative_
($r_{i}>0$) or _obligate_ ($r_{i}<0$). The case of _predation/parasitism_
($b_{ij}<0$, $b_{ji}>0$, i.e. species $j$ benefits at the expense of species
$i$) requires the prey or parasitised species to obtain resources from the
environment ($r_{j}>0$) for the system to be sustainable. When $b_{ij}\approx
0$, this sort of interaction is usually referred to as _commensalism_.
Finally, _competition_ ($b_{ij}<0$, $b_{ji}<0$) occurs when both species are a
hindrance to each other (which of course needs $r_{i}>0$ and $r_{j}>0$ for the
system to be sustainable).

Table 1: Signs of $r_{i}$ and $b_{ij}$ that describe typical ecological interactions. $r_{i}$ | $r_{j}$ | $b_{ij}$ | $b_{ji}$ | _type of interaction_ | abbreviation
---|---|---|---|---|---
$+$ | $+$ | $+$ | $+$ | facultative-facultative mutualism | FFM
$+$ | $-$ | $+$ | $+$ | obligate-facultative mutualism | OFM
$-$ | $-$ | $+$ | $+$ | obligate-obligate mutualism | OOM
$+$ | $+$ | $-$ | $-$ | competition | CMP
$+$ | $+$ | $-$ | $+$ | facultative predator-prey / parasitism | FPP
$+$ | $-$ | $-$ | $+$ | obligate predator-prey / parasitism | OPP
$+$ | any | $\sim 0$ | $+$ | commensalism | COM
Figure 1: Sketch of the model for two microbial species. Figures of microbes
and molecules are purely symbolic, included only for the sake of illustration.
They are meant to represent two species taking different resources from a pool
and interchanging residues. Image created with [24].

Parameters $r_{i}$ and $b_{ij}$ are purely phenomenological. However, in a
real system, they depend on the specific ways in which the species relate with
the environment and with each other. Because of this, the parameters are not
independent. However, to specify how they are connected, we need to delve into
the details of the system we want to model. This is a necessary step to
construct an evolutionary model of an ecosystem [11].

Let us assume that the two microbial species described by Eq. (1) struggle to
survive in a fixed environment (described by the parameters $a_{i}$ and
$c_{i}$), with a given number $K$ of available resources. Each species
consumes $n_{i}$ of these $K$ resources, $q$ of which are common to both
species (shared resources). As a consequence of internal metabolic reactions,
species $i$ produces $m_{i}$ byproducts, $m_{ji}$ of which are useful to the
other species $j$ and $\ell_{i}=m_{i}-m_{ji}$ are not (we do not care about
metabolic waste in this model). Even though $m_{i}$ depends on the details of
the metabolism of species $i$, for the sake of simplicity we can just assume
that it is proportional to the number of consumed resources, i.e.

$m_{1}=\gamma(n_{1}+m_{12}),\qquad m_{2}=\gamma(n_{2}+m_{21}),$ (2)

where $\gamma$ quantifies _metabolic efficiency_. A fraction of species $i$’s
metabolic byproducts will be used by species $j$; hence,

$m_{21}=\alpha_{2}(n_{1}+m_{12}),\qquad m_{12}=\alpha_{1}(n_{2}+m_{21}),$ (3)

where $\alpha_{j}$ ($0<\alpha_{j}<\gamma$) measures _cross-feeding efficiency_
—parameter $\alpha_{i}$ determines the benefit that species $i$ obtains from
species $j$—so it directly influences coefficient $b_{ij}$ in Eq. (1). In this
toy model we assume that these cross-feeding efficiencies can evolve, which
makes the coefficients in Eq. (1) evolve accordingly.

Equation (3) is a linear system whose solution is

$m_{12}=\frac{\alpha_{1}(n_{2}+\alpha_{2}n_{1})}{1-\alpha_{1}\alpha_{2}},\qquad
m_{21}=\frac{\alpha_{2}(n_{1}+\alpha_{1}n_{2})}{1-\alpha_{1}\alpha_{2}},$ (4)

and substituting in Eq. (2),

$m_{1}=\frac{\gamma(n_{1}+\alpha_{1}n_{2})}{1-\alpha_{1}\alpha_{2}},\qquad
m_{2}=\frac{\gamma(n_{2}+\alpha_{2}n_{1})}{1-\alpha_{1}\alpha_{2}}.$ (5)

These equations express the total number of metabolic residues [c.f. Eq. (5)]
and the number of beneficial byproducts [c.f. Eq. (4)] as a function of the
number of resources.

Finally, we need to connect the demographic parameters with this flux of
resources and byproducts. In principle, if resources are available in limited
amounts, microbes would fare better the more resources they have, whereas
metabolic byproducts are costly. It has been argued though that if resources
are abundant becoming a specialist could be advantageous over being a
generalist [25]. The reason is that diversifying mechanisms to use different
resources has a cost in fitness. For the sake of simplicity, we will stick to
the first scenario. Thus, we posit that the growth rate increases proportional
to the amount of consumed resources and decreases with the amount of metabolic
byproducts, i.e.

$\begin{split}r_{1}&=\frac{r_{10}}{K}(n_{1}-m_{1})\\\
&=\frac{r_{10}}{1-\alpha_{1}\alpha_{2}}\big{[}(1-\gamma-\alpha_{1}\alpha_{2})u_{1}-\gamma\alpha_{1}u_{2}\big{]},\\\
r_{2}&=\frac{r_{20}}{K}(n_{2}-m_{2})\\\
&=\frac{r_{20}}{1-\alpha_{1}\alpha_{2}}\big{[}(1-\gamma-\alpha_{1}\alpha_{2})u_{2}-\gamma\alpha_{2}u_{1}\big{]},\end{split}$
(6)

denoting $u_{i}\equiv n_{i}/K$. Notice that the factor $K$ is introduced for
convenience, so as to express everything in terms of the re-scaled variables
$0\leqslant u_{i}\leqslant 1$. (See [26] for a list of the parameters and
symbols used in this model.)

Likewise, the interaction coefficients $b_{ij}$ increase with $m_{ij}$, the
number of byproducts of species $j$ that are useful to species $i$, and
decrease with the competition for the shared resources $q$. Hence,

$\begin{split}b_{12}&=\frac{b_{10}}{K}(m_{12}-q)\\\
&=b_{10}\left[\frac{\alpha_{1}}{1-\alpha_{1}\alpha_{2}}(\alpha_{2}u_{1}+u_{2})-w\right],\\\
b_{21}&=\frac{b_{20}}{K}(m_{21}-q)\\\
&=b_{20}\left[\frac{\alpha_{2}}{1-\alpha_{1}\alpha_{2}}(\alpha_{1}u_{2}+u_{1})-w\right],\end{split}$
(7)

where we have denoted $w\equiv q/K$ ($0\leqslant
w\leqslant\min(u_{1},u_{2})$).

Parameters $r_{i0}$ and $b_{i0}$ are simple dimensional constants that also
set the time scale of the differential equations.

Thus, after these simplifying assumptions, we end up with a model in which the
demographic parameters are expressed in terms of the number of consumed
resources, $n_{1}$, $n_{2}$, the resource sharing $q$, and the cross-feeding
efficiencies $\alpha_{1}$, $\alpha_{2}$. A sketch of the model, where all this
interactions are summarized, is shown in Figure 1. The pool of resources is
represented by the rectangle $K$ in blue, where $n_{1}$ and $n_{2}$ are the
sets of resources useful for species $1$ and $2$, respectively, and $q$ is the
set of shared resources. Residues excreted by each species are in the upper
gray rectangles. Note that the shape and color of microbes are meaningless,
and so are the molecules depicted. They are used just for illustration
purposes.

##  III Adaptive dynamics

According to Eq. (1), the dynamics of species $i$’s population fits the
pattern

$\dot{N}_{i}=N_{i}f_{i}(\Omega,N_{1},N_{2}),$ (8)

where $\Omega$ denotes the set of parameters
$\\{u_{1},u_{2},w,\alpha_{1},\alpha_{2}\\}$. The function
$f_{i}(\Omega,N_{1},N_{2})$ describes the per-capita growth rate, or
_fitness,_ of species $i$, a magnitude that depends on the population growth
parameters as well as the population sizes of the two species involved. Any
steady state of the community, $N_{i}=N_{i}^{*}$, is defined by the equations

$f_{i}(\Omega,N_{1}^{*},N_{2}^{*})=0.$ (9)

When a community is in a given steady state, a mutant may appear in one of the
populations. If the mutant belongs to species $i$, its fitness will be a
function $\bar{f}_{i}(\Omega^{\prime},\Omega,N_{1},N_{2})$, where the prime
parameters are those of the mutant. If this fitness is negative the mutant
will go extinct, otherwise it can increase its frequency in the population and
replace the original genotype.

Adaptive dynamics (AD) is a method to exploit this idea to devise a set of
differential equations for the demographic parameters. Under the assumption
that the parameters of the mutant are a small perturbation of those of the
original genotype, [11] derived these differential equations from the master
equation of the underlying stochastic mutation-selection process.

The derivation of these so-called canonical equations of AD implicitly assumes
that the demographic parameters can vary independently of each other, and that
mutations have the same probability to increase or decrease these parameters
by the same small amount. Neither of these two conditions are met in our model
under the scheme of mutations of this system (increasing or decreasing the
number of resources). One kind of mutations amounts to changing the status of
a randomly chosen external resource, i.e., adding the resource if it is new,
or dropping it if the species was already using it. This means that $n_{i}\to
n_{i}+1$ with a probability proportional to the number of new resources, and
$n_{i}\to n_{i}-1$ with a probability proportional to the number of resources
already in use. But at the same time, $q$, the number of common resources to
both species, can increase, decrease, or remain unchanged. Table 2 summarises
all mutational scenarios and their corresponding probabilities.

Table 2: Set of mutations in the system that change the number of external resources and their corresponding probabilities. Here $n_{i}$ are the number of external resources consumed by species $i$, $q$ the amount of them common to both species, and $u_{i}\equiv n_{i}/K$, $w\equiv q/K$. initial state: $n_{1}$, $q$ | initial state: $n_{2}$, $q$
---|---
mutated state |  | probability | mutated state |  | probability
$n_{1}+1$, | $q+1$ |  | $u_{2}-w$ |  $n_{2}+1$, | $q+1$ |  | $u_{1}-w$
$n_{1}+1$, | $q$ |  | $1-u_{1}-u_{2}+w$ |  $n_{2}+1$, | $q$ |  | $1-u_{1}-u_{2}+w$
$n_{1}-1$, | $q$ |  | $u_{1}-w$ |  $n_{2}-1$, | $q$ |  | $u_{2}-w$
$n_{1}-1$, | $q-1$ |  | $w$ |  $n_{2}-1$, | $q-1$ |  | $w$

The second kind of mutations will change $m_{ij}$, the amount of byproducts
that species $i$ take from species $j$, by $\pm 1$. As $\alpha_{i}$ is a
multiplicative factor,

$\frac{\delta m_{ij}}{m_{ij}}=\frac{\delta\alpha_{i}}{\alpha_{i}},$ (10)

so $\delta m_{ij}=\pm 1$ implies

$\delta\alpha_{i}=\frac{\pm\alpha_{i}}{m_{ij}},$ (11)

or, using Eq. (4),

$\begin{split}\delta\bar{\alpha}_{1}\equiv\frac{u_{2}+\alpha_{2}u_{1}}{1-\alpha_{1}\alpha_{2}}\delta\alpha_{1}=\pm\frac{1}{K},\\\
\delta\bar{\alpha}_{2}\equiv\frac{u_{1}+\alpha_{1}u_{2}}{1-\alpha_{1}\alpha_{2}}\delta\alpha_{2}=\pm\frac{1}{K}.\end{split}$
(12)

On the other hand, $\delta\alpha_{i}<0$ with probability
$m_{ij}/m_{i}=\alpha_{i}/\gamma$ and $\delta\alpha_{i}>0$ with probability
$\ell_{i}/m_{i}=(\gamma-\alpha_{i})/\gamma$. (See [26] for a list of variables
and symbols used in this model .)

In summary, both species can experience two kinds of mutations. The first one
changes the number of resources the species consumes, that is, it may add one
new resource or dispose of an old one, hence increasing or decreasing $n_{i}$
one unit. Notice that this added or removed resource may or not be shared with
the other species. The second kind of mutation is related with its ability to
use byproducts excreted to the environment by the other species—an ability
accounted for by the parameter $\alpha_{i}$.

Reconstructing the procedure of [11] to obtain the canonical equation of AD,
we arrive at a system of differential equations for $u_{1}$, $u_{2}$, $w$,
$\alpha_{1}$, $\alpha_{2}$ (c.f. Eq. (24) and Eq. (28) in Appendix B; see the
Supplementary Material for a more detailed derivation of the equations).
Notice that the stationary populations $N_{1}^{*}$, $N_{2}^{*}$ appear
explicitly in the equations, and that these populations depend in turn on the
parameters $u_{1}$, $u_{2}$, $w$, $\alpha_{1}$, $\alpha_{2}$.

##  IV Numerical simulations

As the differential system (24) is very difficult to discuss analytically, we
have performed numerical simulations of the time evolution of the system for
different initial conditions and sets of parameters. Two timescales are
involved here: the population changes in a fast timescale, whereas the
parameters evolve slowly. The evolutionary parameters are obtained by
integrating Eqs. (24) and (28). We do that using a fourth-order Runge-Kutta
method (RK4) to obtain the parameters
$\\{u_{1}(t),u_{2}(t),w(t),\alpha_{1}(t),\alpha_{2}(t)\\}$ in evolutionary
time. Note that RK4 needs the slope of the function at 4 different times and
that slope depends explicitly on the stationary populations $N_{1}^{*}$,
$N_{2}^{*}$ through Eq.(16), so at each of those times we need to calculate
the stationary populations—which depend, in turn, on the evolutionary
parameters $\\{u_{1},u_{2},w,\alpha_{1},\alpha_{2}\\}$ at those same times.

All simulations have been run with a time-step $\Delta t=10^{-5}$ and were
stopped as soon as they reached a stationary state—when the evolutionary
parameters did not change more than an error tolerance $10^{-4}$ over $10^{4}$
time steps—or after $10^{7}$ times steps.

##  V Evolutionary attractors

Using different sets of parameters and initial conditions, all the
evolutionary attractors that we have observed fit within just a few patterns.
The triad of parameters $(u_{1},u_{2},w)$ (fraction of resources consumed by
each species and fraction of shared resources) is always found to end up in
one of the three forms $(1,u,u)$, $(1-u,u,0)$, or
$\big{(}1-u,u,\min(1-u,u)\big{)}$, with $0\leq u\leq 1$. If $(1,u,u)$ is
reached, the evolutionary attractor turns out to be either competition or
antagonism. In particular, in the case $(1,1,1)$ the attractor is always
competition. As expected, when the final state is mutualism or commensalism no
resources are shared ($w=0$) (there is no need for competing for resources),
so when $(1-u,u,0)$ is reached the system ends up being mutualistic.

With respect to the cross-feeding efficiencies, $\alpha_{i}$, the system
generally evolves towards extreme values of these parameters ($\alpha_{i}=0$
or $\alpha_{i}=\gamma$). Only when the system evolves to competition can some
intermediate values be found. As expected, cross-feeding efficiencies reach
their maximum ($\alpha_{1}=\alpha_{2}=\gamma$) for mutualism and their minimum
($\alpha_{1}=\alpha_{2}=0$) mainly for competition.

It is worth mentioning that all kinds of mutualisms are found to be
evolutionary attractors for some initial conditions. Another relevant
observation is that evolution sometimes drives ecosystems to extinction. This
is no longer a surprise because it is a result that has already been
empirically observed (e.g. in viruses [27]), but the idea that evolution can
degrade ecosystems would have shocked evolutionists of the 19th and early 20th
century because it contradicts the notion of evolution as a sort of
‘optimiser’.

##  VI Evolutionary transitions between types of ecological interactions

In order to illustrate the kind of evolutionary transitions between types of
ecological interactions that this system can produce, we performed an
exhaustive exploration of parameters. We fixed the environmental parameters
$a_{1}=a_{2}=c_{1}=c_{2}=0.001$ and, without loss of generality, chose
$r_{10}=1$ (this simply sets the evolutionary time scale). For the other
species, we explored uniformly the interval $0<r_{20}\leq 1$ and then zoomed
in the region $0<r_{20}\ll r_{10}$ by sampling the interval $0<r_{20}\leq
0.1$. For $b_{i0}$ we sampled uniformly the range $0\leq b_{i0}\leq 0.01$ and
then zoomed in the interval $0\leq b_{i0}\leq 0.001$. Likewise, nine
different, uniformly spaced values of the metabolic efficiency $\gamma$ within
the range $0<\gamma<1$ were explored.

For each set of parameters we generated random initial conditions for $u_{i}$,
$w$, and $\alpha_{i}$ within the ranges $0\leq u_{i}\leq 0.99$, $0\leq
w\leq\min(u_{1},u_{2})$, and $0\leq\alpha_{i}\leq\gamma$, and kept only those
that generated viable populations. Then we let each of these remaining cases
evolve according to Eq. (24) and Eq. (28). We recorded 1000 runs that resulted
in viable populations for each set of parameters, discarding all initial
conditions that led to no viable stationary populations but keeping track of
those that eventually led to extinction. Notice that the number of resources
or the mutation probability are only relevant to set the evolutionary time
scale (see Appendix B).

In order to catalogue the resulting evolutionary attractors we have followed
the classification of Table 1—considering a state as commensalist if one of
the parameters $b_{ij}$ is positive and the other one is smaller than
$10^{-8}$.

The results of this numerical study of the model are summarized in Fig. 2, 3,
4, and 5 (for the distribution of the initial and final states see the
Supplementary Material). In what follows we describe in more detail the
transitions from an initial ecological state to another one that we observed,
depending on the choice of parameters $r_{i0}$ and $b_{i0}$.

###  VI.0.1 Parameters: $r_{20}\leq 1,\ b_{i0}\leq 0.01$

Figure 2: Evolutionary transitions for values of $\gamma$ ranging from $0.1$
to $0.9$. In rows, we show the initial ecological interaction states
(mutualism (MUT), antagonism (ANT), and competition (CMP)), and in columns,
the final ecological interaction states—with the addition of commensalism
(COM) and extinction (EXT). Figures denote the probability of the
corresponding transitions. Parameter values: $r_{20}\leq 1$ and $b_{i0}\leq
0.01$. For $\gamma\leq 0.2$, the initial states are mainly ANT and CMP, and
they mostly end up in MUT or EXT. For $\gamma\geq 0.3$, the initial states are
mainly MUT and remain MUT.

(See Fig. 2.) For this choice of parameters most evolutionary pathways ended
up in mutualism, from $\sim 40\%$ to more than $90\%$ of the cases. For lower
values of $\gamma$, from $0.1\leq\gamma\leq 0.3$, $\sim 40\%$ of cases begin
as antagonism or competition and extinction occurs in more than $20\%$ of all
transitions. For higher values of $\gamma$, from $0.4\leq\gamma\leq 0.7$,
mutualism is the main initial state representing more than $30\%$ of the cases
and extinction reduces to less than $20\%$. For $0.8>\gamma$, mutualism
accounts for $\sim 90\%$ of the initial and final states, and extinction
represent less than $16\%$ of the latter. In most cases, except for
$\gamma=0.9$, commensalism accounts for $\sim 10\%$ of the final states.

###  VI.0.2 Parameters: $r_{20}\leq 0.1,\ b_{i0}\leq 0.01$

Figure 3: As Fig. 2, for parameters values $r_{20}\leq 0.1$ and $b_{i0}\leq
0.01$. For $\gamma\leq 0.2$, the initial states are mainly CMP and ANT,
although they mostly end up in EXT. For $\gamma\geq 0.3$, the initial states
are mainly MUT and remain MUT.

(See Fig. 3.) For this choice of parameters most evolutionary pathways ended
up in mutualism or went extinct. Extinction accounts for more than $45\%$ of
final states for $\gamma\leq 0.3$ reducing its importance for higher values
until $\gamma=0.9$, where represents less than $2\%$. Mutualism begins being
less than $15\%$ of all final states but ended up being almost $90\%$ of the
cases. Commensalism remains representing $\sim 10\%$ of the final states for
all $\gamma$. Initial states are mainly antagonism and competition for lower
values of $\gamma$, representing more than $80\%$ for $\gamma=0.1$ but turned
to mutualism for their higher values, where it represents more than $80\%$ for
$\gamma=0.9$.

###  VI.0.3 Parameters: $r_{20}\leq 1,\ b_{i0}\leq 0.001$

Figure 4: As Fig. 2, for parameters values $r_{20}\leq 1.0$ and $b_{i0}\leq
0.001$. For $\gamma\leq 0.3$, the initial states are mainly CMP, and they
mostly stay as CMP. For $0.3\leq\gamma\leq 0.7$, communities can also begin as
MUT or ANT. For $\gamma\geq 0.6$, most populations end up in CMP or EXT.

(See Fig. 4.) For this parameters the system evolved mainly to competition or
went extinct. With such a small interaction parameters $b_{i0}$ mutualism is
dramatically hindered—even in the cases where one third of the initial states
were mutualistic (for $\gamma=0.5$). More than $40\%$ of the system started in
competition and remained as such or went extinct—the proportion of which
changed upon increasing $\gamma$ from $46.8\%$ competition vs. $24.8\%$
extinction, to $33.3\%$ competition vs. $44.2\%$ extinction.

###  VI.0.4 Parameters: $r_{20}\leq 0.1,\ b_{i0}\leq 0.001$

Figure 5: As Fig. 2, for parameters values $r_{20}\leq 0.1$ and $b_{i0}\leq
0.001$. For $\gamma\leq 0.3$, the initial states are predominantly CMP, and
they mostly end in EXT. For $\gamma\geq 0.4$, the initial states are
predominantly MUT, although EXT continues to be the most common evolutionary
fate. However, for $\gamma\geq 0.8$, COM and ANT become relatively common
ending states.

(See Fig. 5.) As in the previous case, the system cannot evolve into a
mutualistic state, but commensalism accounts for almost $10\%$ of the final
states when $\gamma\leq 0.3$ and more than $15\%$ for $0.4\leq\gamma$.
Antagonism as well accounts for more than $30\%$ of the final cases when
$\gamma=0.9$, being lower for $\gamma\leq 0.8$. Extinction is the main
transition for all cases. Even though competition is the main initial state
for lower values of $\gamma$, being more than $60\%$ of the initial cases for
$\gamma=0.1$ and around $20\%$ for $0.5\leq\gamma$, it shares the distribution
with mutualism and antagonism when $\gamma$ increases, both of them being more
than $20\%$ for all $0.2\leq\gamma$.

##  VII Transient states

Because of the rich evolutionary dynamics of this model, evolutionary
transitions are not the only relevant feature to focus on. Particular
sequences of transient ecological states along evolutionary pathways are as
interesting—if not more so. In Figs. 6, 7, 8, 9, 10, and 11 we show the time
evolution of the demographic parameters of just a few examples, chosen because
of their peculiar sequence of intermediate transitions or because they
resemble actual transitions observed in nature. (For the time evolution of the
stationary populations $N_{i}$ and their evolutionary parameters, see the
Supplementary Material).

Fig. 6 shows two species with positive intrinsic growth rate that are
initially in competition. Over time one of them starts parasitizing the other
until it becomes dependent on it. Eventually, the other species ‘learns’ to
take advantage of the parasite and the relation ends up as a mutualism.
Transitions from antagonistic relations to mutualism are well known to occur
in microbial communities due to environmental pressure and phenotypic
plasticity [22]. However, [28] showed that mutualism arises as an evolutionary
change in a controlled experiment involving Escherichia coli and Salmonella
enterica. E. coli strains went through genetic changes due to a mutation that
led them from generating acetate—a costless byproduct that was useful to
Salmonella—to secreting the more useful galactose, even when this was a costly
byproduct that reduced its intrinsic growth rate. This transition arose when
Salmonella strains were forced and selected to produce methionine, which was
beneficial to E. coli. Even though E. coli mutant populations did not replaced
the ancestral strains—as they were also benefited on this stage—this case
shows how species must pay a cost in order to become mutualists.

Figure 6: Evolutionary transition of the population parameters $r_{i}$ and
$b_{i}$ of two species. The community starts in competition (CMP) and ends in
obligate-facultative mutualism (OFM), going through intermediate states of
facultative parasitism (FPP) and obligate parasitism (OPP). Initial
parameters: $r_{01}=1$, $r_{02}=0.65727$, $b_{01}=0.00414$, $b_{02}=0.00447$,
$u_{1}=0.17967$, $u_{2}=0.78500$, $w=0.10050$, $\alpha_{1}=0.08365$,
$\alpha_{2}=0.32459$, $\gamma=0.5$.

In Fig. 7 species 2 starts as a facultative parasite of species 1. Soon the
parasite becomes dependent, and this situation remains like this for a while,
until suddenly (in evolutionary terms) the relation evolves into a mutualism,
facultative at first for species 1, but eventually obligate for both species.

A similar co-evolutionary pathway can be found in nature between
Macrotermitinae (species 1) and fungi (species 2). According to [29] and [30],
both the fungi and the fungus-growing termites are obligate mutualists since
they need each other in order to survive and reproduce. As stated by [31], it
is plausible that the origin of the termite-fungi mutualistic relation was an
infection (a specialized infestation) of the termites guts with fungi spores,
which led them to defend themselves by domesticating the fungi, controlling
and limiting their growth.

Figure 7: Evolutionary transition of the population parameters $r_{i}$ and
$b_{i}$ of two species from facultative parasitism (FPP) to obligate-obligate
mutualism (OOM), with intermediate states of obligate parasitism (OPP) and
obligate-facultative mutualism (OFM). Initial parameters: $r_{01}=1.0$,
$r_{02}=0.64396$, $b_{01}=0.00474$, $b_{02}=0.00262$, $u_{1}=0.60097$,
$u_{2}=0.62973$, $w=0.16897$, $\alpha_{1}=0.13282$, $\alpha_{2}=0.54215$ and
$\gamma=0.6$.

Fig. 8 is an example of an opposite transition, in which one of the species of
an initially mutualistic system evolves into a parasite of the other (first
facultative, eventually obligate) which, over time, develops a commensalistic
relation with the parasitized species.

Commensalism arises frequently in nature. For instance some algae and
Ascomycota fungi do not form lichens, but descend of lichen-forming fungi
ancestors. Non-lichenized fungi might obtain their nutrients acting as
commensalists of lichen-forming fungi and algae [32]. Although the evolution
of the lichen symbiosis is believed to have appeared multiple independent
times [33], no specific route for this formation has been described with
certainty, to our knowledge.

There are many other cases of known mutualist relations that had become
commensalistic. Zooxanthellae and octocorals form an ancestral facultative
mutualistic relation where the dinoflagellates contribute to the energy budget
of the invertebrates, by being host inside them. However, [34] showed that
there is compelling evidence of a FFM$\rightarrow$COM transition, within the
octocoral family Melithaeidae. Zooxanthellae and antheopleura form a similar
relation and [35] stated that many species of the sea anemones seem to have
lost the mutualistic relationship with the algae. Gymnodinium algae provide an
energetic supply to many marine invertebrates, which in return act as
protective hosts. [36] show, through a molecular phylogenetic analysis, that
Gymnodinium is a genus that incorporate both mutualistic and independent
species even though they all descend from a common symbiotic ancestor.

Similar cases can be found within microbial organisms. In genera Entamoeba and
Trypanosoma, transitions from mutualism to what in our framework can be
described as mutual commensalism have been reported. Clark and Roger [37]
reported that Entamoeba histolytica shows evidence of mitochondrial relics,
which might mean that those organelles, first acquired through endosymbiotic
mechanisms, were eventually lost as a result of an evolutionary process.
Hannaert et al. [38] also reported that two genera of Trypanosomatidae
(Leishmania and Trypanosoma) acquired their plastids probably from mutualistic
algae, and they also lost them leaving only a few residual genes. Such
interpretations are allowing a better understanding of eukaryotic lineages
since organisms formerly classified as divergent ancestors, like some
archezoan protists [37] or nematodes having and lacking their Wolbachia
pipientis simbionts [39], are being understood now as much closer relatives.
These studies show evolutionary transitions that are fundamentally different
those shown in microbial studies due to environmental pressures [22] even
where genetic manipulation is involved [40].

Figure 8: Evolutionary transition of the population parameters $r_{i}$ and
$b_{i}$ of two species from facultative-facultative mutualism (FFM) to
commensalism (COM), with intermediate states of facultative parasitism (FPP)
and obligate parasitism (OPP). Initial parameters: $r_{01}=1.0$,
$r_{02}=0.09827$, $b_{01}=0.00463$, $b_{02}=0.00313$, $u_{1}=0.79063$,
$u_{2}=0.80318$, $w=0.00460$, $\alpha_{1}=0.22106$, $\alpha_{2}=0.12717$, and
$\gamma=0.7$.

In Fig. 9 an initially parasitic relation, with a long period of stasis,
eventually evolves into an interdependent mutualistic relation after crossing
a brief period of facultative mutualism. Notice that the co-evolution of
_Legionella jeonii_ and _Amoeba proteus_ described in the Introduction
illustrates this kind of evolutionary pathway [2].

Figure 9: Evolutionary transition of the population parameters $r_{i}$ and
$b_{i}$ of two species from obligate parasitism (OPP) to obligate-obligate
mutualism (OOM), with an intermediate state of obligate-facultative mutualism
(OFM). The initial OPP has a long stasis, wherear the intermediate OFM is much
shorter in evolutionary scale. Initial parameters: $r_{01}=1$,
$r_{02}=0.04230$, $b_{01}=0.00962$, $b_{02}=0.00096$, $u_{1}=0.45854$,
$u_{2}=0.26628$, $w=0.12217$, $\alpha_{1}=0.28321$, $\alpha_{2}=0.35339$, and
$\gamma=0.6$.

Another example, Fig. 10, exhibits a case that begins and ends in a
mutualistic state, but not before going through a period of parasitism. This
case reveals that, even though some of the cases reported in the statistics of
the previous section appear not to have undergone any transition whatsoever,
they nevertheless can come across different intermediate states before
reaching their evolutionary stable state.

Figure 10: Evolutionary transition of the population parameters $r_{i}$ and
$b_{i}$ of two species from obligate-facultative mutualism (OFM) to obligate-
obligate mutualism (OOM), with an intermediate state of obligate parasitism
(OPP) and obligate-facultative mutualism (OFM). Initial parameters:
$r_{01}=1.0$, $r_{02}=0.80306$, $b_{01}=0.00561$, $b_{02}=0.00310$,
$u_{1}=0.26642$, $u_{2}=0.17242$, $w=0.01485$, $\alpha_{1}=0.40480$,
$\alpha_{2}=0.31492$, and $\gamma=0.6$.

As a last example, Fig. 11 exhibits a case that begins in facultative
parasitism and ends in a obligate-facultative mutualism, going through a
period of facultative-facultative mutualism. This case is similar to the case
of ants and aphids from several genera, whose relations have gone through all
mutualistic and commensalistic nuances, as we pointed out in the Introduction.

Attine ants and fungi might have also undergone a similar pathway, since they
form mutualistic relations that range from almost mutually obligatory to
facultative—at least for the fungi. In particular, [29] indicate that the
fungal symbionts of the higher attines depend almost exclusively on the
ants—although occasionally they might reproduce sexually. However, [41] and
[42] also show that fungi are capable of living without the symbiotic relation
with the ants.

Several hypotheses have been drawn to explain the evolutionary origin of the
attine ant-fungus relation. According to [43], even though it is widely
accepted that fungi were part of the ancestral ant diet (which would
constitute a predatory relation), neutral coexistence, where ants acted as
accidental vehicles of fungal dispersion, may be a more viable explanation (in
which case, the ancestral relation would be commensalistic, instead of
predatory).

Figure 11: Evolutionary transition of the population parameters $r_{i}$ and
$b_{i}$ of two species from facultative parasitism (FPP) to obligate-
facultative mutualism (OFM), with an intermediate state of facultative-
facultative mutualism (FFM). Initial parameters: $r_{01}=1.0$, $r_{02}=1.0$,
$b_{01}=0.00713$, $b_{02}=0.00309$, $u_{1}=0.38958$, $u_{2}=0.96352$,
$w=0.05325$, $\alpha_{1}=0.10848$, $\alpha_{2}=0.08866$, and $\gamma=0.2$.

##  VIII Discussion

Admittedly, the model we have proposed and analyzed in this paper is a far cry
from any realistic description of microbial interactions. For instance, we
have deliberately ignored any detail on the metabolic processes involved—which
are determinant in deciding which products can or cannot be re-used—or kept
the environmental conditions constant—hence neglecting any effect that a
change in the environment might bring into the interactions [21, 44]. For
these and many other drastic simplifications we have assumed in its design, it
is only fair to call our model a ‘toy model’. And yet this is precisely one of
its virtues, because what this model makes clear is that very few assumptions
about the way two species can interact lead to drastic changes in the
ecological scenario. More and more, empirical evidence is showing that
evolutionary transitions in the ecological interactions between species of a
community are the norm, rather than the exception. The model proposed here
provides a proof of concept in this sense, because it links natural
microscopic interactions between microbes to subsequent changes in their
ecological relationships.

It has been hypothesized that environmental fluctuations are one of the forces
that might drive the evolution of microbial consumption and production
capabilities. Our simple model does not dispute this point of view, but the
fact that the random evolution of the interspecific interactions in a constant
environment, subject to some trade-offs, is able to produce such a plethora of
complex transitions between different ecological regimes, can only mean that
the environment may not be their only driver. As a matter of fact,
endosymbiosis is a process that occurs within a constant environment—and one
that our simple model is able to capture. Certainly, environmental changes may
trigger this evolutionary changes, but the subsequent adaptation that they
bring about is due to random changes in the interactions between the species
involved, and may occur even if the environment remains constant.

With a simple bookkeeping of the amount of resources that are taken from the
environment, of those that are excreted as byproducts of metabolic reactions,
and of the amount of byproducts from the other species that can be re-used for
their own purpose, the model can cover virtually all possible ecological
interactions between two species, and show, using adaptive dynamics, that
these interactions evolve, going through different scenarios until reaching a
final stable state. The model also provides some clues about general trends
that can occur in real situations. For instance, there is a marked trend
toward the emergence of mutualistic interactions, even in systems that start
in competition or show antagonistic relationships. Also, many of the observed
pathways have a real counterpart, because similar ones have been documented
for actual species (microbial or otherwise).

The model relies upon the availability of a phenomenological model of
ecological interactions [19, 20] that is capable of describing mutualistic,
competitive, and antagonistic interactions with a simple tuning of the
parameters—in a way that generalized Lotka-Volterra models are not capable of.
This general logistic model of population dynamics is one of the simplest
models capable of exhibiting evolutionary transitions between different
ecological regimes. The clue in devising our present evolutionary model has
been linking those parameters to microscopic interactions between the species
involved—which we have chosen as microbes for the sake of simplicity.

There are at least two ways in which this work can be extended. One is making
more detailed and realistic assumptions on the microscopic interactions that
occur between the two species and see whether the trends observed in this toy
model are or are not confirmed. Obviously, introducing further details will
impose constraints that the present model is currently free of. These
constraints will bias the distribution of scenarios we observe in different
ways, and will do so differently for different species—for which the details
can vary from instance to instance. This will provide interesting information
about the connection between microscopic mechanisms and ecological
transitions.

The other way to extend this work is to consider more than two species. We can
only imagine the richness of ecological scenarios that such an extension will
reveal, even with a mechanism as simple as the straightforward bookkeeping of
resources we have implemented. The computational complexities of this
extension are evident by just looking at the involved calculations that only
two species has led to. This is one of the reasons why we have decided to
postpone such a study for future research—the other one being that just two
species are enough to provide the proof of concept we aimed at with this work.

###### Acknowledgements.

This research has been funded by the Spanish Ministerio de Ciencia, Innovación
y Universidades-FEDER funds of the European Union support, under projects
BASIC (PGC2018-098186-B-I00, J.A.C., and PGC2018-093854-B-I00, J.G. and
J.M.P.), and EVA (CGL2016-77377-R, J.M.I).

##  Appendix A Generalised logistic model of population dynamics

The idea behind the population model of [20] is to extend Velhurst’s logistic
equations of populations

$\dot{N}_{i}=N_{i}\left(\bar{r}_{i}-\bar{a}_{i}N_{i}\right),\quad i=1,\dots
p,$ (13)

by making the parameters $\bar{r}_{i}$ and $\bar{a}_{i}$ to depend on the
interactions with the environment as well as the population sizes of all
species in the community as

$\bar{r}_{i}=r_{i}+\sum_{k=1}^{p}b_{ik}N_{k},\qquad\bar{a}_{i}=a_{i}+c_{i}\sum_{k=1}^{p}b_{ik}N_{k}.$
(14)

Here $r_{i}$ is the intrinsic growth rate of species $i$, $b_{ik}$ is the rate
of benefit (if positive) or hindrance (if negative) on species $i$ due to the
interaction with species $k$, and $p$ is the total number of species in the
ecosystem. The coefficients $a_{i}$ measure intraspecific competitions (hence
$a_{i}>0$) due to a limitation of the environmental resources. As a matter of
fact, in the standard Velhurst’s model $a_{i}^{-1}$ directly measures the
carrying capacity of the environment. As of $c_{i}$, the effect of these
coefficients is better seen if we rewrite Eq. (13) and Eq. (14) as

$\dot{N}_{i}=N_{i}\left(r_{i}-a_{i}N_{i}\right)+N_{i}\left(1-c_{i}N_{i}\right)\sum_{k=1}^{p}b_{ik}N_{k}.$
(15)

We can clearly see in this expression that choosing $c_{i}>0$ induces a
saturation on the interaction of the community with the focal species.

For only two species the stationary solutions can be obtained as:

$\begin{split}N^{\ast}_{1}=\frac{-a_{1}a_{2}+b_{12}b_{21}+b_{21}c_{2}r_{1}-b_{12}c_{1}r_{2}\pm\sqrt{\Delta}}{2b_{21}(b_{12}c_{1}+a_{1}c_{2})},\\\
N^{\ast}_{2}=\frac{-a_{1}a_{2}+b_{12}b_{21}-b_{21}c_{2}r_{1}+b_{12}c_{1}r_{2}\pm\sqrt{\Delta}}{2b_{12}(b_{21}c_{2}+a_{2}c_{1})},\end{split}$
(16)

where

$\begin{split}\Delta=&\,4b_{21}(b_{12}c_{1}+a_{1}c_{2})(a_{2}r_{1}+b_{12}r_{2})\\\
&+[a_{1}a_{2}-b_{21}(b_{12}+c_{2}r_{1})+b_{12}c_{1}r_{2}]^{2}.\end{split}$
(17)

According to [20], the linear stability of the finite stationary solution for
two species ($N^{\ast}_{1},N^{\ast}_{2}$) can be analyzed from the Jacobian
matrix

$J(N^{\ast}_{1},N^{\ast}_{2})=\left(\begin{array}[]{cc}-r_{1}-b_{12}N^{\ast}_{2}&b_{12}N^{\ast}_{1}\left[1-c_{1}N^{\ast}_{1}\right]\\\
b_{21}N^{\ast}_{2}\left[1-c_{2}N^{\ast}_{2}\right]&-r_{2}-b_{21}N^{\ast}_{1}\end{array}\right)$
(18)

##  Appendix B Adaptive dynamics equations for the two-species ecosystem

According to Eq. (1), the per-capita fitness of species $i$ is given by

$f_{i}(\Omega,N_{1},N_{2})=r_{i}(\Omega)-a_{i}N_{1}+b_{ij}(\Omega)(1-c_{i}N_{i})N_{j},$
(19)

with $r_{i}(\Omega)$ and $b_{ij}(\Omega)$ given by Eq. (6) and Eq. (7)
respectively. If the parameters of the mutant are $u^{\prime}_{i}=u_{i}+\delta
u_{i}$, $w^{\prime}=w+\delta w$, with $\delta u_{i}=\pm 1/K$ and $\delta
w=0,\pm 1/K$ (see Table 2), then—given the linearity of $f_{i}$ with respect
to these parameters—the per-capita fitness of the mutant in a steady state
community can be written as

$\bar{f}_{i}(u^{\prime}_{i},w^{\prime},\Omega,N_{1}^{*},N_{2}^{*})=f_{i,u}\delta
u_{i}-f_{i,w}\delta w,$ (20)

where

$\begin{split}f_{i,u}&=\frac{r_{i0}(1-\gamma-\alpha_{1}\alpha_{2})+b_{i0}\alpha_{1}\alpha_{2}(1-c_{i}N_{i}^{*})N_{j}^{*}}{1-\alpha_{1}\alpha_{2}},\\\
f_{i,w}&=b_{i0}(1-c_{i}N_{i}^{*})N_{j}^{*}.\end{split}$ (21)

Accordingly, the canonical equations of AD for this sort of mutants read (see
[11])

$\begin{split}\dot{u}_{1}&=\mu N_{1}^{*}\sum_{\delta u_{1},\delta w}\delta
u_{1}p_{1}(\delta u_{1},\delta w)\big{[}\bar{f}_{1}(u_{1}+\delta
u_{1},w+\delta w,\Omega,N_{1}^{*},N_{2}^{*})\big{]}_{+},\\\ \dot{u}_{2}&=\mu
N_{2}^{*}\sum_{\delta u_{2},\delta w}\delta u_{2}p_{2}(\delta u_{2},\delta
w)\big{[}\bar{f}_{2}(u_{2}+\delta u_{2},w+\delta
w,\Omega,N_{1}^{*},N_{2}^{*})\big{]}_{+},\\\ \dot{w}&=\mu
N_{1}^{*}\sum_{\delta u_{1},\delta w}\delta wp_{1}(\delta u_{1},\delta
w)\big{[}\bar{f}_{1}(u_{1}+\delta u_{1},w+\delta
w,\Omega,N_{1}^{*},N_{2}^{*})\big{]}_{+}\\\ &\phantom{=}+\mu
N_{2}^{*}\sum_{\delta u_{2},\delta w}\delta wp_{2}(\delta u_{2},\delta
w)\big{[}\bar{f}_{2}(u_{2}+\delta u_{2},w+\delta
w,\Omega,N_{1}^{*},N_{2}^{*})\big{]}_{+},\end{split}$ (22)

where the sums run over all the corresponding mutations and $p_{i}(\delta
u_{i},\delta w)$ are their respective probabilities, according to Table 2, and
$\mu$ is the probability of mutation per reproduction event. The function
$[x]_{+}$ stands for $x$ if $x\geq 0$ and $0$ otherwise.

Substituting Eq. (20) and Eq. (21), performing the sums, and re-scaling
evolutionary time with $2K^{2}/\mu$ yields the set of differential equations

$\begin{split}\dot{u}_{1}&=2N_{1}^{*}\left\\{[f_{1,u}-f_{1,w}]_{+}(u_{2}-w)-w[f_{1,w}-f_{1,u}]_{+}\right.\\\
&\phantom{=}\left.+[f_{1,u}]_{+}(1-u_{1}-u_{2}+w)-[-f_{1,u}]_{+}(u_{1}-w)\right\\},\\\
\dot{u}_{2}&=2N_{2}^{*}\left\\{[f_{2,u}-f_{2,w}]_{+}(u_{1}-w)-w[f_{2,w}-f_{2,u}]_{+}\right.\\\
&\phantom{=}\left.+[f_{2,u}]_{+}(1-u_{1}-u_{2}+w)-[-f_{2,u}]_{+}(u_{2}-w)\right\\},\\\
\dot{w}&=2N_{1}^{*}\left\\{[f_{1,u}-f_{1,w}]_{+}(u_{2}-w)-w[f_{1,w}-f_{1,u}]_{+}\right\\}\\\
&\phantom{=}+2N_{2}^{*}\left\\{[f_{2,u}-f_{2,w}]_{+}(u_{1}-w)-w[f_{2,w}-f_{2,u}]_{+}\right\\},\end{split}$
(23)

which, using the identity $[x]_{+}=(x+|x|)/2$, can be rewritten as

$\begin{split}\dot{u}_{1}&=N_{1}^{*}\left\\{f_{1,u}-f_{1,w}u_{2}+|f_{1,u}-f_{1,w}|(u_{2}-2w)\right.\\\
&\phantom{=}\left.+|f_{1,u}|(1-u_{2}-2u_{1}+2w)\right\\},\\\
\dot{u}_{2}&=N_{2}^{*}\left\\{f_{2,u}-f_{2,w}u_{1}+|f_{2,u}-f_{2,w}|(u_{1}-2w)\right.\\\
&\phantom{=}\left.+|f_{2,u}|(1-u_{1}-2u_{2}+2w)\right\\},\\\
\dot{w}&=N_{1}^{*}\left\\{(f_{1,u}-f_{1,w})u_{2}+|f_{1,u}-f_{1,w}|(u_{2}-2w)\right\\}\\\
&\phantom{=}+N_{2}^{*}\left\\{(f_{2,u}-f_{2,w})u_{1}+|f_{2,u}-f_{2,w}|(u_{1}-2w)\right\\}.\end{split}$
(24)

On its side, mutants that change their parameter $\alpha_{i}$ have a per-
capita fitness

$\bar{f}_{i}(\alpha^{\prime}_{i},\Omega,N_{1}^{*},N_{2}^{*})=f_{i,\alpha}\delta\bar{\alpha}_{i}+O\left(\delta\bar{\alpha}_{i}^{2}\right),$
(25)

where $\delta\bar{\alpha}_{i}=\pm 1/K$ and

$f_{i,\alpha}=\frac{b_{i0}(1-c_{i}N_{i}^{*})N_{j}^{*}-r_{i0}\gamma}{1-\alpha_{1}\alpha_{2}}.$
(26)

Accordingly,

$\begin{split}\dot{\alpha}_{1}&=\mu
N_{1}^{*}\sum_{\delta\alpha_{1}}\delta\alpha_{1}p_{1}(\delta\alpha_{1})\big{[}\bar{f}_{1}(\alpha_{1}+\delta\alpha_{1},\Omega,N_{1}^{*},N_{2}^{*})\big{]}_{+},\\\
\dot{\alpha}_{2}&=\mu
N_{2}^{*}\sum_{\delta\alpha_{2}}\delta\alpha_{2}p_{2}(\delta\alpha_{2})\big{[}\bar{f}_{2}(\alpha_{2}+\delta\alpha_{2},\Omega,N_{1}^{*},N_{2}^{*})\big{]}_{+},\end{split}$
(27)

which, neglecting $O\left(\delta\bar{\alpha}_{i}^{2}\right)$ terms and re-
scaling time as before, becomes

$\begin{split}\dot{\alpha}_{1}&=\frac{2N_{1}^{*}}{\gamma}\left\\{[f_{1,\alpha}]_{+}(\gamma-\alpha_{1})-[-f_{1,\alpha}]_{+}\alpha_{1}\right\\}\\\
&=N_{1}^{*}\left\\{f_{1,\alpha}+|f_{1,\alpha}|\frac{\gamma-2\alpha_{1}}{\gamma}\right\\},\\\
\dot{\alpha}_{2}&=\frac{2N_{2}^{*}}{\gamma}\left\\{[f_{2,\alpha}]_{+}(\gamma-\alpha_{2})-[-f_{2,\alpha}]_{+}\alpha_{2}\right\\}\\\
&=N_{2}^{*}\left\\{f_{2,\alpha}+|f_{2,\alpha}|\frac{\gamma-2\alpha_{2}}{\gamma}\right\\}.\end{split}$
(28)

## References

  * Park  _et al._ [2004] M. Park, S. Yun, M. Kim, J. Chun, and T. Ahn, Phylogenetic characterization of Legionella-like endosymbiotic X-bacteria in Amoeba proteus: a proposal for ‘Candidatus Legionella jeonii’ sp. nov., Environ. Microbiol. 6, 1252 (2004).
  * Jeon and Lorch [1967] K. W. Jeon and I. J. Lorch, Unusual intra-cellular bacterial infection in large, free-living amoebae, Exp. Cell Res. 48, 236 (1967).
  * Jeon [1992] K. W. Jeon, Macromolecules involved in the amoeba-bacteria symbiosis, J. Protozool. 39, 199 (1992).
  * Jeon [1995] K. W. Jeon, Bacterial endosymbiosis in amoebae, Trends Cell. Biol. 5, 137 (1995).
  * Machado  _et al._ [2001] C. A. Machado, E. Jousselin, F. Kjellberg, S. G. Compton, and E. A. Herre, Phylogenetic relationships, historical biogeography and character evolution of fig-pollinating wasps, Proc. R. Soc. B 268, 685 (2001).
  * Shingleton and Stern [2002] A. W. Shingleton and D. L. Stern, Molecular phylogenetic evidence for multiple gains or losses of ant mutualism within the aphid genus Chaitophorus, Mol. Phylogenet. Evol. 26, 26 (2002).
  * Sakata [1994] H. Sakata, How an ant decides to prey on or to attend aphids, Res. Popul. Ecol. 36, 45 (1994).
  * Stadler and Dixon [2005] B. Stadler and A. F. G. Dixon, Ecology and evolution of aphid-ant interactions, Annu. Rev. Ecol. Evol. Syst. 36, 345 (2005).
  * Offenberg [2001] J. Offenberg, Balancing between mutualism and exploitation: the symbiotic interaction between Lasius ants and aphids, Behav. Ecol. Sociobiol. 49, 304 (2001).
  * Sachs and Simms [2006] J. L. Sachs and E. L. Simms, Pathways to mutualism breakdown, Trends Ecol. Evol. 21, 585 (2006).
  * Dieckmann and Law [1996] U. Dieckmann and R. Law, The dynamical theory of coevolution: a derivation from stochastic ecological processes, J. Math. Biol. 34, 579 (1996).
  * Dercole and Rinaldi [2008] F. Dercole and S. Rinaldi,  _Analysis of Evolutionary Processes: The Adaptive Dynamics Approach and Its Applications_ (Princeton University Press, Princeton, New Jersey, 2008).
  * Doebeli [2011] M. Doebeli,  _Adaptive Diversification_ (Princeton University Press, Princeton, New Jersey, 2011).
  * Lotka [1925] A. J. Lotka,  _Elements of Physical Biology_ (Williams and Wilkins Company, Baltimore, 1925).
  * Volterra [1926] V. Volterra, Fluctuations in the abundance of a species considered mathematically, Nature 118, 558 (1926).
  * Volterra [1928] V. Volterra, Variations and Fluctuations of the Number of Individuals in Animal Species Living Together, ICES J. Mar. Sci. 3, 3 (1928).
  * May [1981] R. M. May, Models for two interacting populations, in  _Theoretical Ecology: Principles and Applications_ , edited by R. M. May and A. R. McLean (Oxford Univesity Press, Oxford, UK, 1981) pp. 78–104.
  * Wright [1989] D. H. Wright, A simple, stable model of mutualism incorporating handling time, Am. Nat. 134, 664 (1989).
  * García-Algarra  _et al._ [2014] J. García-Algarra, J. Galeano, J. M. Pastor, J. M. Iriondo, and J. J. Ramasco, Rethinking the logistic approach for population dynamics of mutualistic interactions, J. Theor. Biol. 363, 332 (2014).
  * Stucchi  _et al._ [2020] L. Stucchi, J. M. Pastor, J. García-Algarra, and J. Galeano, A general model of population dynamics accounting for multiple kinds of interaction, Complexity 2020, 10.1155/2020/7961327 (2020).
  * Thompson [1988] J. N. Thompson, Variation in interspecific interactions, Annu. Rev. Ecol. Syst. 19, 65 (1988).
  * Hoek  _et al._ [2016] T. A. Hoek, K. Axelrod, T. Biancalani, E. A. Yurtsev, J. Liu, and J. Gore, Resource availability modulates the cooperative and competitive nature of a microbial cross-feeding mutualism, PLOS Biology 15, e1002606 (2016).
  * Tilman [2004] D. Tilman, Niche tradeoffs, neutrality, and community structure: A stochastic theory of resource competition, invasion, and community assembly, Proc. Natl. Acad. Sci. (USA) 101, 10854 (2004).
  * BioRender [2019] BioRender, (2019), https://biorender.com, last accessed on 2019-07-10.
  * Schwartz and Hoeksema [1998] M. W. Schwartz and J. D. Hoeksema, Specialization and resource trade: biological markets as a model of mutualisms, Ecology 79 (1998).
  * [26] See Supplemental Material at [URL to be inserted by publisher] for details about the derivation of the model as well as supplementary figures.
  * Turner and Chao [1999] P. E. Turner and L. Chao, Prisoner’s dilemma in an RNA virus, Nature 398, 441 (1999).
  * Harcombe  _et al._ [2018] W. R. Harcombe, J. M. Chacón, E. M. Adamowicz, L. m. Chubiz, and C. J. Marx, Evolution of bidirectional costly mutualism from byproduct consumption, PNAS 115, 12000 (2018).
  * Nobre  _et al._ [2010] T. Nobre, C. Rouland-Lefevre, and D. K. Aanen, Biology of termites: a modern synthesis (Springer, Dordrecht, 2010) Chap. Comparative Biology of Fungus Cultivation in Termites and Ants, pp. 193–210.
  * Aanen  _et al._ [2002] D. K. Aanen, P. Eggleton, C. Rouland-Lefevre, T. Guldberg-Frøslev, S. Rosendahl, and J. J. Boomsma, The evolution of fungus-growing termites and their mutualistic fungal symbionts, PNAS 99, 14887 (2002).
  * Margulis and Sagan [2002] L. Margulis and D. Sagan,  _Acquiring Genomes: A Theory of the Origin of Species_ (Basic Books, New York, NY, 2002).
  * Lutzoni  _et al._ [2001] F. Lutzoni, M. Pagel, and V. Reeb, Major fungal lineages are derived from lichen symbiotic ancestors, Nature 411, 937 (2001).
  * Divakar  _et al._ [2015] P. K. Divakar, A. Crespo, M. Wedin, S. D. Leavitt, D. L. Hawksworth, L. Myllys, B. McCune, T. Randlane, J. W. Bjerke, Y. Ohmura, I. Schmitt, C. G. Boluda, D. Alors, B. Roca-Valiente, R. Del-Prado, C. Ruibal, K. Buaruang, J. Núñez-Zapata, G. A. de Paz, V. J. Rico, M. C. Molina, J. A. Elix, T. L. Esslinger, I. K. K. Tronstad, H. Lindgren, D. Ertz, C. Gueidan, L. Saag, K. Mark, G. Singh, F. D. Grande, S. Parnmen, A. Beck, M. N. Benatti, D. Blanchon, M. Candan, P. Clerc, T. Goward, M. Grube, B. P. Hodkinson, J.-S. Hur, G. Kantvilas, P. M. Kirika, J. Lendemer, J.-E. Mattsson, M. I. Messuti, J. Miadlikowska, M. Nelsen, J. I. Ohlson, S. Pérez-Ortega, A. Saag, H. J. M. Sipman, M. Sohrabi, A. Thell, G. Thor, C. Truong, R. Yahr, D. K. Upreti, P. Cubas, and H. T. Lumbsch, Evolution of complex symbiotic relationships in a morphologically derived family of lichen-forming fungi, New Phytol. 208, 1217 (2015).
  * Oppen  _et al._ [2005] M. J. H. V. Oppen, J. . C. . Mieog, C. A. Sánchez, and K. E. Fabricius, Diversity of algal endosymbionts (zooxanthellae) in octocorals: the roles of geography and host relationships, Mol Ecol 14, 2403 (2005).
  * Geller and Walton [2001] J. B. Geller and E. D. Walton, Breaking up and getting together: evolution of symbiosis and cloning by fission in sea anemones (genus Anthopleura), Evolution 55, 1781 (2001).
  * Wilcox [1998] T. P. Wilcox, Large-Subunit ribosomal RNA systematics of symbiotic dinoflagellates: morphology does not recapitulate phylogeny, Mol. Phylogenet. Evol. 10, 436 (1998).
  * Clark and Roger [1995] C. G. Clark and A. J. Roger, Direct evidence for secondary loss of mitochondria in Entamoeba histolytica, PNAS 92, 6518 (1995).
  * Hannaert  _et al._ [2003] V. Hannaert, E. Saavedra, F. Duffieux, J.-P. Szikora, D. J. Rigden, P. A. M. Michels, and F. R. Opperdoes, Plant-like traits associated with metabolism of Trypanosoma parasites, PNAS 100, 1067 (2003).
  * Casiraghi  _et al._ [2004] M. Casiraghi, O. Bain, R. Guerrero, C. Martin, V. Pocacqua, S. L. Gardner, A. Franceschi, and C. Bandi, Mapping the presence of Wolbachia pipientis on the phylogeny of filarial nematodes: evidence for symbiont loss during evolution, Int. J. Parasitol. 34, 191 (2004).
  * LaSarre  _et al._ [2017] B. LaSarre, A. L. McCully, J. T. Lennon, and J. B. McKinlay, Microbial mutualism dynamics governed by dose-dependent toxicity of cross-fed nutrients, ISME Journal 11, 337 (2017).
  * Schultz and Brady [2008] T. R. Schultz and S. G. Brady, Major evolutionary transitions in ant agriculture, PNAS 105, 5435 (2008).
  * Currie  _et al._ [2003] C. R. Currie, B. Wong, A. E. Stuart, T. R. Schultz, S. A. Rehner, U. G. Mueller, G.-H. Sung, J. W. Spatafora, and N. A. Straus, Ancient tripartite coevolution in the attine ant-microbe symbiosis, Science 299, 386 (2003).
  * Mueller  _et al._ [2001] U. G. Mueller, T. R. Schultz, C. R. Currie, R. M. M. Adams, and D. Malloch, The origin of the attine ant-fungus mutualism, Q Rev Biol 76, 169 (2001).
  * Hernández [1998] M. J. Hernández, Dynamics of transitions between population interactions: A nonlinear interaction $\alpha$-function defined, P. Roy. Soc. B-Biol. Sci. 265, 1433 (1998).