Coupled fast and slow feedbacks lead to continual evolution: A general modeling approach

TL;DR

This paper presents a general modeling framework demonstrating that the interplay of fast positive and slow negative feedbacks can lead to continual evolution, emphasizing robustness and minimal assumptions in evolutionary dynamics.

Contribution

It introduces a minimal, general model showing how coupled feedbacks induce continual evolution, advancing understanding of biotic interaction-driven evolutionary processes.

Findings

Coupled feedbacks cause continual evolution under certain timescale separations.

The model demonstrates robustness of continual evolution to perturbations.

Minimal conditions are sufficient for ongoing evolutionary change.

Abstract

The Red Queen Hypothesis, which suggests that continual evolution can result from solely biotic interactions, has been studied in macroevolutionary and microevolutionary contexts. While the latter has been effective in describing examples in which evolution does not cease, describing which properties lead to continual evolution or to stasis remains a major challenge. In many contexts it is unclear which assumptions are necessary for continual evolution, and whether described behavior is robust under perturbations. Our aim here is to prove continual evolution under minimal conditions and in a general framework, thus automatically obtaining robustness. We show that the combination of a fast positive and a slow negative feedback causes continual evolution with a single evolving trait, provided the ecological timescale is sufficiently separated from the timescales of mutations and negative feedback. Our approach and results form a next step towards a deeper understanding of the evolutionary dynamics resulting from biotic interactions.

Equations

$$\frac{du}{dt}$$

$$\frac{d\phi }{dt}$$

$$\frac{dA}{dt}$$

$$\frac{dB}{dt}$$

$$\frac{d\phi }{dt}$$

$$\frac{d{u}_i}{dt}$$

$$\frac{d\phi }{dt}$$

$$R=\frac{A}{A+B},$$

$$\frac{dR}{dt}$$

$$\frac{d\phi }{dt}$$

$$\frac{\partial f}{\partial R}>0.$$

$$R=\frac{A}{A+B}.$$

$$\frac{dA}{dt}$$

$$\frac{dB}{dt}$$

$$\frac{d\phi }{dt}$$

$$\left(\frac{\partial }{\partial A}-\frac{\partial }{\partial B}\right)f_A>0,\mathrm{and}\left(\frac{\partial }{\partial B}-\frac{\partial }{\partial A}\right)f_B>0.$$

$$f_A(A,B,\phi )>0$$

$$f_A(A,B,\phi )<0$$

$$f_B(A,B,\phi )>0$$

$$f_B(A,B,\phi )<0$$

$$\frac{\partial f_A}{\partial A}(A_0,B,\phi )\le 0,$$

$$\frac{\partial f_A}{\partial B}<\frac{\partial f_A}{\partial A}.$$

$$\frac{\partial {\eta }_A}{\partial B}=-\frac{\frac{\partial f_A}{\partial B}}{\frac{\partial f_A}{\partial A}}$$

$$\left(\frac{\partial }{\partial A}-\frac{\partial }{\partial B}\right)h>0.$$

$$f_B>f_A$$

$$f_A>f_B$$

$$A\cdot f_A=0\mathrm{and}B\cdot f_B=0.$$

$$\frac{dA}{dt}$$

$$\frac{dB}{dt}$$

$$\frac{dA}{dt}=\frac{dB}{dt}=0,$$

$$M:=\frac{{\sum }_{i=0}^{1}i{u}_i}{{\sum }_{i=0}^1{u}_i}.$$

$$\frac{d{u}_i}{dt}=u_i\left({\mu }_i\left(1-\frac{{\sum }_{i=0}^1{u}_i}{k}\right)-d\right),$$

$$\frac{dM}{dt}$$

$$=\frac{\left({\sum }_{i=0}^{1}i\frac{d{u}_i}{dt}\right)\left({\sum }_{i=0}^1{u}_i\right)-\left({\sum }_{i=0}^{1}i{u}_i\right)\left({\sum }_{i=0}^1\frac{d{u}_i}{dt}\right)}{{\left({\sum }_{i=0}^1{u}_i\right)}^2}.$$

$$\frac{{\sum }_{i=0}^{1}i\frac{d{u}_i}{dt}}{{\sum }_{i=0}^{1}i{u}_i}>\frac{{\sum }_{i=0}^1\frac{d{u}_i}{dt}}{{\sum }_{i=0}^1{u}_i}.$$

Full text

1. Introduction

Species typically live in complex ecosystems with many interactions among them and external drivers. The evolutionary dynamics of a species in a complex ecosystem can be caused by the properties of the species, by the interaction with the coexisting species or by external drivers. To better understand to what extent the emerging ecological and evolutionary dynamics is caused by intrinsic properties of the species or biotic interactions (within- and between species interaction) we need to ignore abiotic drivers. Such an abiotically unchanging environment may lead to a static adaptive landscape, where the adaptation will follow a path towards a peak in that landscape—reachable or not [1]. Since the major part of any individual’s environment typically is composed of other (evolving) species, any species’ environment will in general change even without external abiotic variation, both through ecological and evolutionary changes. Hence, the adaptive landscape will be dynamic, which can lead to continual co-evolutionary dynamics, where traits over evolutionary time are fluctuating, escalating or chasing each other [2].

Of particular interest to us is the emergence of the so-called Red Queen Dynamics [3], a concept that has had a major influence on micro- and macroevolutionary theory. [4] and others [5, 6, 7] have analyzed a macro-evolutionary model aiming at understanding under what conditions continual Red Queen dynamics and stasis results from within-system biotic interactions: Stenseth and Maynard Smith demonstrated that both could result depending upon the nature of the within system biotic interactions — without being able translated into ecological terms what these conditions were. [8] found that symmetric competitive interactions are more likely to lead to stasis. Another approach is studying interactions of a few species in more mechanistic detail, microevolutionary Red Queen (RQ) dynamics [2].

Theoretical studies of micro-evolutionary RQ dynamics mostly use methods based on adaptive dynamics and quantitative genetics [9, 10, 11, 12, 13, 14]. The adaptive dynamics approach assumes the ecological dynamics have reached an equilibrium and studies the invasion of individuals with a slightly deviating phenotype of the adaptive trait. The advantage of this approach is that it allows for a rigid theoretical analysis of the system. The quantitative genetics approach does not assume the ecological dynamics to be in equilibrium, but has a timescale separation between the evolutionary adaptation and the ecological dynamics. If the evolutionary rate of change is very slow, the quantitative genetics approach becomes similar to the adaptive dynamics approach. Both methods assume that adaptive traits evolve along a fitness gradient. Studies focussing on predator-prey and host-parasite systems have been able to reach conclusions about conditions that increase or decrease the chance of RQ dynamics in a specific setting (e.g. fast adaptation is less likely to lead to RQ dynamics [10] and RQ dynamics requires an intermediate harvesting efficiency of the prey [9]).

Most of the above mentioned studies use specific functional forms for their analysis, hampering the generalisation of the obtained results. Meta-analysis (such as the one by [15]) can provide some more general insights, but the conclusions are still limited, especially since many studies use similar equations. To obtain general results and therefore a broad understanding of what ecological interactions can cause certain evolutionary patterns, such as continual evolution, we need as general models as possible. With this contribution we aim at extending the theoretical understanding of under which conditions continual evolution and stasis will result. For this purpose we use a very general model with few assumptions regarding the form of the model functions. With such a general approach (making a minimum of assumption on the functional forms and the parameters) we increase the robustness of the obtained results. We find that a system with slow and fast feedback interactions exhibits continual Red Queen type of dynamics depending on the timescales. Moreover we allow for a polymorphic population, not constraining the distribution of phenotypes that may be present in the population, and mutations of small and large effects.

2. Model description and the emerging eco-evolutionary dynamics

2.1. Model description

In order to focus our argument, we use a general model of an evolving trait for a single species, which can be extended to multiple species:

[TABLE]

Here $u$ represents the population density distribution over the trait space. In the general case the evolving trait is not specified and therefore its trait space could be anything. The external factors (abiotic factors, species without evolving traits) are captured in $\varphi$. The function $f(u,\varphi)$ describes the growth of the population (which can depend on the distribution in trait space $u$ as well as the external factors $\varphi$). The function $g(u)$ describes the change in the trait due to mutations. The function $g$ does not incorporate how well these phenotypes perform (the change due to mutations is modeled in an unbiased manor), the change of the trait to more fit individuals comes from the modelling of the growth of all phenotypes. The function $h(u,\varphi)$ describes the change in external factors, which can depend on the phenotype densities $u$. Time scale differences between growth, external factors and mutations are captured in $\epsilon_{e}$ and $\epsilon_{m}$.

We include a possible polymorphic trait distribution (as is preferable argued by [16] and also used in [17]). A polymorphic trait distribution arises easily with asexual reproduction or traits that are determined by a few loci, but can also arise when assortive mating develops (see the discussion in [18]). We use an approach including genetic variation by modeling the populations with differential equations, which also releases the assumption that mutations have to cause an infinitesimal change in the trait value.

An example of a slow negative feedback is a predator prey interaction, where the generation time of the predator is much slower than that of the prey and the predator predates preferable on prey with a specific phenotype. An example of the specificity of the predator is predators preferring a certain size of prey (e.g. bears preferring large salmon [19] or zooplankton preferring phytoplankton with a certain nitrogen:phosphorus ratio [14]). Other examples besides predator-prey interactions include: a herbivore that changes the vegetation which then becomes less suitable for the herbivore; micro-organisms that produce a compound that inhibits them but is diluted in a large volume (to achieve the feedback being slow); a species with a complex life cycle with a habitat shift (with an evolvable component), if a crowded habitat has a delayed effect on the habitat quality; or a more complicated response of the whole ecosystem that leads to a negative feedback (e.g. with intransitive cycles as in [20]). The last example includes cases where the feedback is caused by human intervention, such as flu vaccinations—common viruses will be vaccinated against during the next season.

The fast positive feedback is an Allee effect on the phenotype. This could be due to for example finding mates (if the trait has effect on suitability of mates), cooperation in defense or in feeding. An example is when the adaptive trait is foraging on a certain vegetation type, the positive feedback could be that the more individuals forage in the same place, the better they are protected against predation. Or in the case of a habitat shift in the life cycle, similar habitat shift will share the second habitat with more individuals, providing mates or protection.

2.2. A single species with two phenotypes shows continual evolutionary dynamics when the negative feedback is sufficiently slow

The best way to understand the system is to reduce it to the simplest version. Therefore we start with an simple version of a trait with only two possible phenotypes, one that is interacting with a external feedback and one that is not (or less strongly) (see Figure 1A). Although the main objective of this simplification is to understand the mechanisms that lead to continual evolution, this representation also reflects a biological system where a trait is either present or not at all, e.g. choosing one host or another. When the trait also affects mate choice (e.g. because they will be located near the same host), the trait has a positive feedback. The slow negative feedback results from the fact that one of the hosts can develop defenses.

We analyse the system where both phenotypes and the feedback change over time. We assume that the phenotypes have a positive feedback effect, i.e. the more individuals with a certain phenotype, the higher the growth rate of that phenotype. To be precise, we study the following system with the interacting phenotype $A$ (the phenotype choosing the host that can develop defenses), the (less) interacting phenotype $B$ (the phenotype choosing the host that cannot develop defenses) and the feedback compound $\varphi$ (the amount of defense in the host population that can develop defenses in the example):

[TABLE]

The growth functions of the phenotypes $f_{A}$ and $f_{B}$ include the positive feedback. $\varphi$ is increased more by $A$ than by $B$ which is incorporated in the function $h$. Mutations are possible from one phenotype to the other with the mutation function $g$. $\epsilon_{m}$ and $\epsilon_{e}$ denote the different timescales of the population dynamics (no $\epsilon$), the feedback dynamics ($\epsilon_{e}$) and the evolutionary dynamics ($\epsilon_{m}$). We have modeled a positive effect of $A$ on $\varphi$ and a negative effect of $\varphi$ on $A$, but all our results are also valid for the opposite, then we can substitute $\varphi$ by $-\varphi$.

We prove that, for the system in Figure 1A and Equations (2) and with only slight restrictions on the functions $f$, $g$ and $h$, if we can choose $\epsilon_{e}$ and $\epsilon_{m}$ sufficiently small, the system will exhibit continual evolution for almost all initial states. The detailed proof is in the Supplementary Information section 4.1 and 4.2; the reasoning is as follows: First, we prove the case of a constant population size, where we can describe our system by the feedback compound $\varphi$ and the ratio $(R)$ between the two phenotypes. Observing that the proof for this case does not rely upon knowing the exact value of $\frac{dR}{dt}$ but merely the sign, we can extend the proof to systems where the population size is not fixed, by adopting an assumption controlling the sign of $\frac{dR}{dt}$. There are several possibilities but we used the ’unique stable value assumption’: when we keep $\varphi$ fixed and know $A$ there is only one value of $B$ that guarantees a stable state, and likewise when we know $B$ there is only one value of $A$. The proof only requires knowledge of the sign in a subregion of the domain.

The idea of the proof for a variable population size is illustrated in Figure 2B. The diagrams show an approximation of the dynamics in the ($A$,$B$) phase planes. This approximation is without the mutation term, but since mutations are rare (in our case the mutation term is small), the actual phase plane is very similar to the one shown. When at some point the feedback level is intermediate, there are 4 intersections of the lines $\frac{dA}{dt}=0$ and $\frac{dB}{dt}=0$ of which one is repelling (the origin), one is a saddle point (the point in the middle) and the other two are attracting (the nature of the equilibria can be seen from the sign of the derivatives). Depending on the initial densities of $A$ and $B$ the system goes to an attracting point near the $A$-axis or the $B$-axis. When this is near the $A$-axis, $B$ will be almost 0 and since $A$ is interacting strongly with the feedback compound, $\varphi$ will increase. This will in turn change the phase plane diagram and the intersection near the $A$-axis will become repelling. Then there is only one attracting intersection and that is near the $B$-axis and the system will approach that state. Since the dynamics in $\varphi$ are slower than in $A$ and $B$, the system will come close to the intersection near the $B$-axis before $\varphi$ changes significantly. When $\varphi$ changes we are back to the initial figure for intermediate $\varphi$, but now $\varphi$ will continue to decrease and the intersection near the $B$-axis will disappear. This will continue indefinitely and therefore will lead to continual evolution.

We only need some restrictions on the functions $f$, $g$ and $h$ and do not use any functional forms in our proof. Figure 2 illustrates the simulation of an example system. We show continual evolutionary cycles for 2 different initial states if $\epsilon_{e}$ and $\epsilon_{m}$ are chosen sufficiently small. With a high mutation rate and/or a fast feedback a stable equilibrium will be reached (leading to stasis), while lack of mutations causes one of the two phenotypes to go extinct.

2.3. A multi-valued trait

For studying the case of a multi-valued trait, the evolving trait is represented as a range of values, and the external factor is represented by a scalar quantity (normalised between 0 and 1). We can then describe the system with the following set of differential equations:

[TABLE]

We split the population in discrete groups $u_{i}$ with similar trait values. We show that the assumptions of a fast positive feedback and a slow negative feedback on the phenotype in this otherwise general model lead to continual evolution. These assumptions are given by restrictions on the functions $f$ and $h$. We can extend the reasoning from the previous section to a multi-valued trait. With a multi-valued trait we mean a trait with more than two phenotypes or an approximation of a continuous trait, such as length or level of toxin production. An overview of the system is shown in Figure 3A, where the trait values are scaled from 0 to 1 according to the strength of their negative feedback. We call these phenotypes $i$, and the density of individuals with phenotype $i$ is denoted by $u_{i}$. The feedback compound is again called $\varphi$ and can be the same quantities as in the previous section. For the positive feedback we use the population average of the trait, which we call $M$. We assume that when the level of $\varphi$ is so low that phenotypes with higher $i$ have higher growth rates for a certain phenotype distribution, growth rates will keep increasing with $i$ when we increase $M$ (and vice versa, if $\varphi$ is high and growth rates decrease with $i$, a decrease in $M$ will not change this). We use the assumption that the growth rate is either strictly increasing or strictly decreasing in the trait (therefore one of the extreme traits always has the highest fitness).

Analogous to the system with two phenotypes, we can prove the existence of evolutionary fluctuations under conditions mimicking our earlier assumptions. In the multiple phenotype setting there is one subtle point in the proof, where an additional assumption is needed to guarantee relatively quick changes in phenotypes. One possible additional assumption guaranteeing such rapid fluctuations is the existence of mutations from any trait value to any other trait value. In practical settings this assumptions will not always be applicable. For example, the assumption is more realistic when the trait value describes enzyme levels or levels of toxin production, than when considering length of a species, as we expect length to evolve gradually. An alternative assumption will be discussed later.

The proof that we can choose the timescales to obtain continual evolutionary dynamics is outlined in Fig. 3B, and described in detail in the Supplementary Information section 4.3. Here we provide an outline: We scale $\varphi$ also between 0 and 1 and without loss of generality may assume that $\varphi$ tends to $M$. Then, when for example $M$ is high and $\varphi$ is low, $\varphi$ will increase (blue arrows in Fig. 3B) because $\varphi$ goes to $M$. At some point $\varphi$ and $M$ will cross, in which case $M$ must already be declining and will continue to decline. The assumed fast positive feedback ensures that $M$ decreases drastically in a time interval in which $\varphi$ hardly changes, leading to symmetric conditions and hence continual evolutionary oscillations (green arrows in Fig. 3B).

Again our proof is context free and applies to a wide variety of equations; an example is given in Figure 4. In Figure 4A demonstrates that turning off evolution leads to extinction of all but one phenotype, fast feedback leads to an equilibrium, and slow feedback leads to continual fluctuations. Figure 4B shows the distribution of phenotypes over time, and in 4C it is shown that the diversity of phenotypes increases at shifts of the mean trait value.

If we do not want to assume the existence of mutations from any phenotype to any other phenotype, we can adopt an alternative assumption, namely that the negative feedback is chosen sufficiently slow depending on the rate of mutations. In some circumstances this may even mean that the rate at which $\varphi$ changes is even slower than the mutation rate. This alternative assumption again may or may not be desirable in practice.

2.4. Extensions to multiple feedbacks

In natural systems, multiple traits determine the fitness in an individual. An extension of our model to more traits will give a more realistic picture of the dynamics we expect in natural populations. We extended the system to a species with two different traits that are both associated with a different negative feedback. The number of possible phenotypes is now the product of the number of phenotypes for each of the traits. In the most simple case, there are four phenotypes in the population: (1) high susceptibility for both feedbacks, (2) low susceptibility for both feedbacks, and (3 and 4) high susceptibility for one and low for the other (Figure 5A). Simulations with this extended system show that this system can also lead to continual adaptation and that the dynamics generally become irregular (Figure 5B). Looking at the change of the diversity index over time, we can conclude that there is no fixed period in the dynamics (Figure 5C). The example with multiple traits shows that we do not expect a natural system to come back to exactly the same state in a regular fashion.

3. Discussion

We have shown that a simple motif of fast positive feedback and slow negative feedback leads to continual evolution. Instead of working with a specific model we have shown that continual evolution occurs for other similar models. We allow for phenotypic variation in the population and large effect mutations. We automatically show robustness to model form, because the evolutionary dynamics do not depend on specific equations or parameter values, and also robustness to small fluctuations. Large fluctuations might affect the dynamics, but because our dynamics are reachable from almost all initial conditions, the system will return to continual evolution after a large fluctuation. Previously, it has been shown that continual evolution can be found in specific simple models. In the literature robustness is indicated by showing that a range of parameter values or a different equation will lead to the same result [21, 10]. Here we try to avoid the use of specific equations and parameter values to obtain more general results and stronger robustness. Previous examples of papers not using specific functional forms consider the possibility of the occurrence of evolutionary sliding [21], the conclusion that symmetric interactions are more likely to lead to stasis [8] and the Red Queen evolution in bacterial communities with cyclic inhibition in three species [20], although all of these papers include some restrictions on the models. [21] focusses only on slow-fast models with separated timescales, [8] only consider bilinear species interactions and [20] only consider a fixed set of three species that form a non-transitive cycle with their species interactions.

Our results are accompanied by specific examples, but the results are not restricted to those examples. In all our examples the continual evolutionary dynamics stem from switching between multiple ecological attractors that are steady states (ecogenetically driven Red Queen Dynamics mode B[22]). Our conclusions remain the same when these attractors are limit cycles instead of steady states. An example of this type of dynamics is given in [22] Figure 4 and, interestingly, directly using our assumptions on this example removes the evolutionary oscillations (only ecological cycles remain, with ecologically driven Red Queen Dynamics—that is fluctuations in traits are fast and small and only follow the ecological population density dynamics). However, when we make the negative feedback slower by decreasing the parameters for the predator dynamics tenfold (parameters r4 and $\gamma$), we retrieve the ecogenetically driven Red Queen Dynamics (results not shown).

Some of our results are in line with and extend previous conclusions. The dynamics of the continual evolution are very similar to microevolutionary Red Queen dynamics; the difference is that continual evolution does not require co-evolution. [22] mention that Red Queen dynamics with a single evolving trait are possible when the dynamics are ecologically driven (the traits follow ecological dynamics) or ecogenetically driven and switching between two different ecological attractors. Continual evolution of a single trait might be quite prevalent, since examples include predator-prey systems in which the prey evolves but the predator does not (or much slower; as might be relatively common, see e.g. [23]). Moreover, our two trait example reproduces the result that fast adaptation is less likely to lead to RQ dynamics [10]. A new observation is that for a multi-valued trait, evolution should not be too slow relative to the negative feedback, and that ecology and evolutionary dynamics can influence each other. The rationale behind this is that the average population trait has to change drastically while the level of the negative feedback variable changes relatively little. If a trait evolved gradually that means that evolution should not be too much slower than the dynamics of the feedback, because then the feedback would ’catch up’ with the average population trait while this average is near the equilibrium value. An experimental example of such kind of dynamics is that in cases such as a predator-prey system with a slow predator, the non-optimal predator phenotypes might stay in the population and not go extinct. In that case the species does not need to get a de novo mutations, but can get the new phenotype through selection and recombination. An example is a laboratory system of an algae and a rotifer [24]. This type of RQ dynamics where the traits stay in the population will be more common in higher organisms and is an example of evolutionary rescue. Faster evolutionary change through mutations might also occur when species are regularly exposed to different environments and evolve adaptability [25, 26].

Phenotypic variation is often left out of Red Queen dynamics analysis. Outright restricting phenotypic variation can have a pronounced effect on the results. We have adapted a competitor-competitor model (Figure 6A) with originally an adaptive dynamics approach that previously showed Red Queen Dynamics to an instance with our assumptions (the possibility of a polymorphic trait distribution and larger mutations). As shown in Figures 6B, switching to the modelling approach used in this paper removes the RQ dynamics that were found in the original paper. However, when we change the timescales and make the feedback (competitor species in this example) slower, we retrieve the feedback, in line with the results of this paper (see Figure 6C). In [10], the RQ dynamics remain with phenotypic variation (data not shown), but there the RQ dynamics are ecologically driven: the adaptation follows the ecological predator-prey cycles. Other methods that do include phenotypic varions are the model by [9], which does contain a polymorphic model which, different from our model, is also stochastic; models of evolutionary branching and extinction using adaptive dynamics (summarized in [18]); and [17], where a Lotka-Volterra predator-prey model with a polymorphic population similar to our methods shows examples of RQ dynamics.

Collecting evidence of RQ dynamics from natural systems is a difficult task [15] and in experiments the long term of the measurements is a problem. One experiment that shows the effect of evolution on predator-prey dynamics [10] suggests RQ dynamics, but due to the measurement time it cannot be excluded that the cycles will dampen. While there is an increase in long term adaptation experiments these are usually in bacterial systems while most of the modelling has focused on sexually reproducing predator-prey communities. Here we have tried to include these systems by allowing for phenotypic variation and larger mutation effect sizes. In bacterial systems polymorphic populations are common, even under reasonably constant conditions as recently shown in the Long Term Evolution Experiment with E. coli [27]. The polymorphic trait distribution might result from density dependent dynamics, as shown theoretically in a chemostat [28]. Therefore our results are an important addition for linking theoretical to experimental observations of evolutionary dynamics.

Although we have tried to keep the modelling in this paper general, our results have some limitations. We show that with some assumptions we can guarantee that continual dynamics will arise with enough time scale separation. Systems that do not follow our assumptions (e.g. fast negative feedbacks, no positive feedback) are more likely to lead to stasis. However we do not prove any conditions that necessarily lead to stasis. Although we tried to keep our equations general, we did not include individual variation within the population. It would be interesting to see if we can combine the methods used in this paper with the polymorphic population and link of the ecological and evolutionary timescales with individual-based models. We limited our results here to only one evolving population, but the feedback we describe could come from another evolving species or be a result of a complex ecosystem. A next step would be to extend this mathematical approach to more than one evolving species. The results in this paper are constrained to micro-evolutionary dynamics, it would be interesting to extend these results to macro-evolutionary phenomena. In this line, there is a nice analogy with the results of [29]. They show that most coevolutionary dynamics are found with intermediate diversity, where not all niches are filled. Here we see that if we get diversification, and therefore more different phenotypes, the co-evolutionary dynamics cease.

Without constraining both the functional forms used in our model and the phenotypic diversity within the system we have demonstrated that a fast positive feedback combined with a slow negative feedback always leads to continual dynamics with the proper timescales. By so doing we have improved the understanding of continual evolution and co-evolution in a large class of models, and may be used to predict evolutionary dynamics without knowledge of the exact equations describing the system.

4. Supplementary information to:

Coupled fast and slow feedbacks lead to continual evolution: A general modeling approach

4.1. One species with two phenotypes and a fixed population size

4.1.1. The system

We first consider the simplest possible system where the combination of a slow negative and a fast but weaker positive feedback cause continual evolution. We consider a single species with two strains $A$ and $B$, and we assume that the total population size $A+B$ remains constant. We will write

[TABLE]

hence $R$ takes on values in the interval $[0,1]$. We introduce a second variable $\varphi$ playing the role of a negative feedback, and will be assumed to also take on values in the interval $[0,1]$. We assume that the rates at which the concentrations $A$, $B$ and $\varphi$ vary depends only on these three variables. Since $A+B$ is assumed to remain constant we can reduce the system of differential equations to the following suggestive form

[TABLE]

Here $f,g$, and $h$ are assumed to be continuously differentiable, and the constants $\epsilon_{m}$ and $\epsilon_{e}$ are assumed to be sufficiently small. The function $f$ represents the ecological dynamics of the populations $A$ and $B$. Since an extinct population cannot reproduce we naturally obtain the extinction assumption

1. (EXT)

The function $f$ takes on the value [math] when $R\in\{0,1\}$.

The function $g$ represents mutations from population $A$ to population $B$ and vice versa. While $\varphi$ will influence the ecology, we assume that $\varphi$ does not affect the mutation rates. Since it is natural to assume that the number of mutations increases with the size of the population, we obtain the mutation assumption

1. (MUT)

The function $g$ is decreasing in $R$, $g(0)>0$ and $g(1)<1$.

The constants $\epsilon_{m}$ and $\epsilon_{e}$ represent the different time scales of the dynamical system. By choosing $\epsilon_{m}$ and $\epsilon_{e}$ sufficiently small, we guaranteed that slower time scales for the mutations and for the ecology of $\varphi$ than for the ecology of $A$ and $B$.

We make three further assumptions: a condition on $f$ that guarantees a positive feedback in the $R$-variable, and two conditions that combined guaranteed a stronger negative feedback on $R$ caused by the slowly adapting variable $\varphi$.

1. (PF)

The function $f$ is strictly decreasing in $\varphi$, and

[TABLE] 2. (NF1)

The function $h$ is decreasing in $\varphi$ and increasing in $R$. 3. (NF2)

There exist $0<\varphi_{--}<\varphi_{++}<1$ such that for $\varphi<\varphi_{--}$, the value of $f$ is strictly positive for all $R\in(0,1)$, and for $\varphi>\varphi_{++}$, the value of $f$ is strictly negative for all $R\in(0,1)$. Similarly, given any $\varphi\in(0,1)$ the function $h(R,\varphi)$ is stricly positive for $R$ suficiently large, and strictly negative for $R$ sufficiently small.

We will prove that this combination of assumptions guarantees large fluctuations in the values of $R$ and $\varphi$.

Let $K\subset(0,1)\times(0,1)$ be a closed subset of states $(R,\varphi)$.

Lemma 1**.**

When $\epsilon_{m}$ and $\epsilon_{e}$ are chosen sufficiently small, the orbit of almost every initial condition $(R_{0},\varphi_{0})\in K$ will eventually leave $K$.

Proof.

By assumption (iii), the differential equation will have strictly positive divergence on $K$ by choosing $\epsilon_{m}$ and $\epsilon_{e}$ sufficiently small. In other words, the dynamical system will be strictly area expanding on $K\cap\{R\geq 1\}$ with respect to the standard Euclidean area. It follows that it is not possible for the orbits of a set of positive area to remain in $K$. ∎

Note that for this result to hold on a given set $K$, assumptions (iii) is only necessary on $K$. In fact, in our simulations below these two assumptions will not be globally satisfied.

Lemma 2**.**

By choosing $\epsilon_{m}$ and $\epsilon_{e}$ sufficiently small, almost every orbit converges to an invariant closed curve not completely contained in $K$.

Proof.

By the Poincaré-Bendixson Theorem, every orbit must converge to either a fixed point or a invariant closed curve.

Independently of the values of $\epsilon_{m}$ and $\epsilon_{e}$, we can construct a large rectangle $[a,b]\times[c,d]\subset(0,1)\times(0,1)$ that must contain all fixed points. Indeed, by choosing $c<\varphi_{--}$ and $d>\varphi_{+}+$ it follows that there are no fixed points for $\varphi<c$ or $\varphi>d$. By assumption (NF2) we can choose $a$ sufficiently small and $b$ sufficiently large so that for $\varphi\in[c,d]$ and $R<a$ the variable $\varphi$ must decrease, while for $\varphi\in[c,d]$ and $R>b$ the variable $\varphi$ must increase. Thus, there are no fixed points outside of the rectangle $[a,b]\times[c,d]$.

As remarked above, the dynamical system can be made area expanding in the rectangle by choosing $\epsilon_{m}$ and $\epsilon_{e}$ sufficiently small. Thus, the fixed points in $[a,b]\times[c,d]$ cannot be stable, nor can the rectangle contain an invariant closed curve. ∎

Theorem 3**.**

For $\epsilon_{e}$ and $\epsilon_{m}$ sufficiently small almost every orbit will converge to a periodic orbit that intersects both $R>b$ and $R<a$.

Proof.

Let $0<R_{-}<R_{+}<1$ be such that $h$ is strictly positive for $R>R_{+}$ and $\varphi<\varphi_{++}$, and strictly negative for $R<R_{-}$ and $\varphi>\varphi_{--}$. Since $f$ is strictly increasing in $R$ and strictly decreasing in $\varphi$, the curve $\{f=0\}$ is a strictly increasing graph that does not intersect $\{\varphi>\varphi_{++}\}$ or $\{\varphi<\varphi_{--}\}$. Denote its graph by $\Gamma$, and write $(R_{+},\varphi_{+})$ for the intersection point of $\Gamma$ with $\{R=R_{+}\}$. Let $R_{++}>R_{+}$.

By choosing $\epsilon_{m}$ sufficiently small we can guarantee that $\frac{dR}{dt}>0$ for $R=R_{++}$ and $\varphi<\hat{\varphi}_{+}$, where $\hat{\varphi}_{+}\in(\varphi_{+},\varphi_{++})$. It follows that whenever an orbit intersects the rectangle $[R_{++},1]\times[0,\hat{\varphi}_{+}]$, the orbit will remain in this region for finite time until it crosses the upper boundary in some point $(R_{t},\hat{\varphi}_{+})$. Note that mutations prevent either population from going extinct, so $R_{t}<1$.

At time $t$, the variable $\varphi$ will continue to increase for as long as $\varphi<\varphi_{++}$ and $R\geq R_{+}$. But since $f$ is negative for $\varphi\geq\varphi_{++}$, it follows that the orbit must cross the line $\{R=r_{+}\}$. The first time this occurs must therefore be in a point $(R_{+},\varphi_{t^{\prime}})$ for which $\varphi_{t^{\prime}}>\hat{\varphi}_{+}$.

Since $\hat{\varphi}_{+}$ is strictly larger than $\varphi_{+}$, it follows that $f<0$ for all $R\leq R_{+}$ and all $\varphi>\varphi_{+}^{\prime}$, where $\varphi_{+}^{\prime}\in(\varphi_{+},\hat{\varphi}_{+})$ can be chosen independent of the exact intersection point $(R_{+},\hat{\varphi}_{+})$. By choosing $\epsilon_{e}$ and $\epsilon_{m}$ sufficiently small, the orbit can be guaranteed to reach any sublevel $R<R_{--}<R_{-}$ while $\varphi$ remains at least $\varphi_{+}^{\prime}$. Since our assumptions are completely symmetric, a symmetric argument shows that $R$ remains smaller than $R_{--}$ until $\varphi<\varphi_{-}$, where $(R_{-},\varphi_{-})\in\Gamma$, and later returns to the region $R>R_{++}$.

As remarked earlier, we can make sure that all fixed points are contained in the rectangle $K=[R_{--},R_{++}],[\varphi_{-}^{\prime},\varphi_{+}^{\prime}]$, where $\varphi_{-}^{\prime}<\varphi_{-}$ is defined analogously as $\varphi_{+}^{\prime}$. Since $\varphi_{+}^{\prime}>\varphi_{+}$ and $\varphi_{-}^{\prime}<\varphi_{-}$ can be chosen independently of the points $(R_{+},\hat{\varphi}_{+})$ and $(R_{-},d_{-})$, it follows that the entire forward orbit avoids the rectangle $K$. Thus the Carleson-Bendixson Theorem implies that the orbit must converge to a periodic cycle. By Lemma 1 we can guarantee that almost every orbit in $K$ must eventually leave $K$, which completes the proof. ∎

4.2. One species with two phenotypes and a variable population size

In the previous section we assumed that $\frac{dR}{dt}$ depends only on $R$ and $\varphi$. When the total population size is not assumed to be fixed, such an assumption is not realistic. However, it turns out that the mechanism causing the drastic fluctuations can still occur without this assumption. The important observation is that to prove the existence of large fluctuations, we merely need to know the signs of $\frac{dR}{dt}$ and $\frac{d\varphi}{dt}$, and not the exact values.

Suppose now that we have two distinct traits, with population sizes $A$ and $B$. We do not assume that $A+B$ is constant, but will still write

[TABLE]

Thus $R$ takes on values in the interval $[0,1]$. We assume that the total population will remain bounded, from above as well as below, so that $R\rightarrow 0$ corresponds to the extinction of population $A$, and $R\rightarrow 1$ to the extinction of population $B$. One naturally assumes that apart from mutations $dA/dt=0$ when $A=0$, and similarly $dB/dt=0$ when $B=0$, to obtain a system of differential equations of the form

[TABLE]

As before, $\varphi$ will play the role of a strong negative feedback. We will assume that the functions $f_{A},f_{B},g,$ and $h$ are all continuously differentiable, and the constants $\epsilon_{m}$ and $\epsilon_{e}$ will be chosen arbitrarily small.

The existence of a positive ecological feedback on the populations $A$ and $B$ leads to the following assumption:

1. (PF)

For all values of $A,B,$ and $\varphi$ we have

[TABLE]

Thus, if an amount of $B$ is replaced by an equal amount of $A$ while $\varphi$ remains fixed, the fitness of the population $A$ increases while the fitness of the population $B$ decreases.

The next assumption, which we will refer to as the unique stable value assumption, will be used to draw conclusions about the sign of $\frac{dR}{dt}$ without knowing the exact values of $A$ and $B$. We do not claim that this assumption is necessary, but it turns out to be satisfied in many models and is convenient:

1. (USV1)

For fixed values of $B$ and $\varphi$, there is a unique value $A_{0}\geq 0$ such that

[TABLE]

for $0\leq A<A_{0}$, and

[TABLE]

for $A>A_{0}$. 2. (USV2)

For fixed values of $A$ and $\varphi$, there is a unique value $B_{0}\geq 0$ such that

[TABLE]

for $B<B_{0}$, and

[TABLE]

for $B>B_{0}$.

For $\varphi$ fixed we can therefore consider $A_{0}$ as a graph $\eta_{A}(B,\varphi)$, and similarly $B_{0}=\eta_{B}(A,\varphi)$. The positive feedback condition guarantees the following:

Lemma 4**.**

The derivatives $\frac{\partial\eta_{A}}{\partial B}$ and $\frac{\partial\eta_{B}}{\partial A}$ are both strictly less than $-1$, and possibly $-\infty$.

Proof.

We prove the first statement, the second is analogous. Let $A_{0},B,\varphi$ be a triple for which $f_{A}(A_{0},B,\varphi)=0$. The unique stable value assumption implies that

[TABLE]

and the positive feedback assumption gives

[TABLE]

The function $\eta_{A}$ is implicitly defined by $f_{A}(\eta_{A}(B,\varphi),B,\varphi)=0$. When $\frac{\partial f_{A}}{\partial A}\neq 0$ the Implicit Function Theorem gives

[TABLE]

which gives $\frac{\partial\eta_{A}}{\partial B}<-1$. When $\frac{\partial f_{A}}{\partial A}=0$ we have $\frac{\partial f_{A}}{\partial B}<0$ by assumption (PF), and hence $\frac{\partial\eta_{A}}{\partial B}=-\infty$. ∎

We assume the existence of a strong negative feedback:

1. (NF1)

The difference $f_{A}-f_{B}$ decreases with $\varphi$, while

[TABLE] 2. (NF2)

There exists $0<\varphi_{--}<\varphi_{++}<1$ such that

[TABLE]

whenever $\varphi>\varphi_{++}$, while

[TABLE]

whenever $\varphi<\varphi_{--}$.

Hence, for $\varphi$ sufficiently large and $\varphi$ sufficiently small the level sets $\{f_{A}=0\}$ and $\{f_{B}=0\}$ do not intersect. We may write $\varphi_{++}$ and $\varphi_{--}$ for the maximal resp. minimal value of $\varphi$ for which the level sets intersect. For our argument it will be necessary that both $\varphi_{++}$ and $\varphi_{--}$ can be reached, hence we make a “transitivity” assumption:

1. (Tr)

The function $h$ is strictly positive when $\varphi\leq\varphi_{++}$, $B=0$ and $A=\eta_{A}(0,\varphi_{++})$. Similarly, the function $h$ is strictly negative when $\varphi\geq\varphi_{--}$, $A=0$ and $B=\eta_{B}(0,\varphi_{--})$.

Let us consider, for a fixed value of $\varphi$, joint solutions of the two equations

[TABLE]

It is clear that, besides the origin, there always is a unique solution on each of the axes. For $A,B>0$ it follows from Lemma 4 that there either is no solution or a unique solution, depending on the value of $\varphi$. For $\varphi=\varphi_{++}$ the intersection point of $\{f_{A}=0\}$ and $\{f_{B}=0\}$ lies in the axis $\{A=0\}$, while for $\varphi=\varphi_{--}$ the intersection point lies in the axis $\{B=0\}$. See Figure 8 for a simple depiction of the level sets $\{f_{A}=0\}$ (in red) and $\{f_{B}=0\}$ (in blue).

Let us consider the ecological dynamics in the $(A,B)$-plane caused by the differential equations

[TABLE]

for a fixed value of $\varphi$. When $\varphi\leq\varphi_{--}$ or $\varphi\geq\varphi_{++}$ there are three fixed points where

[TABLE]

namely the origin and the two intersection points of the curves $\{f_{A}=0\}$ and $\{f_{B}=0\}$ with the respective axes $\{B=0\}$ and $\{A=0\}$. The origin is always repelling. One of the points on the axes is a saddle fixed point, with stable manifold equal to the axis. The third fixed point is attracting, and all orbits of initial values not lying on the axes converge to this attracting fixed point.

When $\varphi_{--}<\varphi<\varphi_{++}$ there are four fixed points. Again the origin is a repelling fixed point. There are again two fixed points on the axes, which are now both attracting. Finally, there is an intersection point of the curves $\{f_{A}=0\}$ and $\{f_{B}=0\}$, and assumption (i) implies that this is a hyperbolic saddle fixed point. Its stable manifold is the separatrix of the two attracting basins.

Let us now consider the effect of the mutations, represented by $\epsilon_{m}\cdot g(A,B)$, on the dynamics in the $(A,B)$-plane in the case $\varphi_{--}<\varphi<\varphi_{++}$. The behavior near each of the fixed points is stable under small $C^{1}$-perturbations, and the qualitative behavior of the system is robust. Hence, by choosing $\epsilon_{m}$ sufficiently small, there will still be a repelling fixed point at the origin, and a saddle point with separatrix near the intersection point of the curves $\{f_{A}=0\}$ and $\{f_{B}=0\}$. The rest of the quadrant is attracted to neighborhoods of the original attracting fixed points. These neighborhoods can be chosen arbitrarily small by choosing $\epsilon_{m}$ sufficiently small. Note that the addition of mutations causes the axes to be repelling, hence the attracting fixed point no longer lies on the axis but sufficiently nearby.

Finally let us consider the full three-dimensional dynamical system, taking into account that $\varphi$ is not fixed. We assume that for given $t_{0}$ we have $\varphi(t_{0})<\varphi_{++}$, and that $(A(t_{0}),B(t_{0}))$ lies in the small attracting region near the $A$-axis. For any $\delta>0$ we can, by taking $\epsilon_{e}$ and $\epsilon_{m}$ sufficiently small, assume that $A(t_{0})-\eta_{A}(0,\varphi_{++})>-\delta$ and $B<\delta$. It follows that $h(A(t_{0}),B(t_{0}),\varphi(t_{0}))$ is strictly positive. It follows that $\varphi$ will continue to grow while $\varphi<\varphi_{++}+\delta$ and while $(A(t),B(t))$ remains trapped in the small attracting neighborhood near the $A$-axis.

Let us consider the first time $t_{1}$ at which either $\varphi\geq\varphi_{++}+\delta$ or at which $(A,B)$ leaves the small attracting neighborhood. In either case it follows that the orbit $(A(t),B(t))$ is guaranteed to approach the small attracting neighborhood near the $B$-axis. If $\epsilon_{e}$ is sufficiently small the value of $\varphi$ can only decrease arbitrarily little while this happens. The conclusion is that we end up with a time $t_{2}>t_{1}$ when $\varphi(t_{1})>\varphi_{--}$ and $(A(t_{2}),B(t_{2}))$ lies in the attracting neighborhood near the $B$-axis. By the symmetry of our assumptions the process will repeat itself. We have proved the following.

Theorem 5**.**

Let $0<a<b<1$, $\varphi_{--}<c<d<\varphi_{++}$, and write $K$ for the rectangle $[a,b]\times[c,d]$. Then for $\epsilon_{e}$ and $\epsilon_{m}$ sufficiently small there exists an orbit $(A(t),B(t),\varphi(t)$ for which the coordinates $(R(t),\varphi(t))$ avoid $K$, and for which $R(t)$ fluctuates between values larger than $b$ and smaller than $a$.

4.3. One species with a trait with multiple phenotypes

Let us recall that in the first setting, where we considered two traits and a constant total population size, the argument used that the rate at which the ratio $R=\frac{A}{A+B}$ changes is determined only by $R$ and $\varphi$. When the assumption on the constant population size was dropped, this was no longer the case, but we saw that the main idea of the argument can still be used because in the appropriate regions in $(R,\varphi$)-coordinates the sign of $\frac{\partial R}{\partial t}$ is still known. One can easily imagine that, under appropriate assumptions, the same approach can be used when dealing with more than two phenotypes. The good choice of $R$, a projection from many variables onto $\mathbb{R}$, is of course essential to the argument.

Let us present a particular family of systems with finitely many trait values. We consider $n+1$ different strains $u_{i}$, where the index $i\in[0,1]$ is of the form $\frac{j}{n}$. As before, we will stipulate a fast positive feedback and a slow but dominating negative feedback caused by a variable $\varphi$. We will assume that $\varphi$ is produced more by strains $u_{i}$ for $i$ large, and for fixed values of $u=(u_{0},\ldots,u_{1})$ converges towards the average index $M$ given by

[TABLE]

In the setting with multiple traits, knowledge of $M$ and $\varphi$ will not be sufficient to determine $dM/dt$ and $d\varphi/dt$. However, we will present assumptions where the signs of $dM/dt$ and $d\varphi/dt$ can be determined, at least in the relevant regions, which will be sufficient to guarantee large fluctuations.

Looking first only at the ecological factors of the dynamical system, ignoring the effects due to mutations, we assume that the densities of the individuals with phenotype $i$ change according to the equations

[TABLE]

where the growth rates $\mu=(\mu_{0},\ldots,\mu_{1})$ are strictly positive, depend continuously on $u$ and $\varphi$, and for any given values of $u$ and $\varphi$ the tuple $(\mu_{0},\cdots,\mu_{1})$ is assumed to be either independent of $i$ or strictly monotonic in $i$. When $\mu$ is constant the relative population sizes do not change, and it follows that $\frac{dM}{dt}=0$.

Lemma 6**.**

When $\mu$ is strictly increasing (resp. decreasing) in $i$, the rate $\frac{dM}{dt}$ is strictly positive (resp. negative).

Proof.

Let us assume that $\mu$ is strictly increasing, the argument is identical when $\mu$ is decreasing. We note that

[TABLE]

Since we are only interested in the sign of $\frac{dM}{dt}$, we can drop the denominator and hence need to prove that

[TABLE]

Plugging in the formula for $\frac{du_{i}}{dt}$, we note that the terms $(1-\frac{\sum_{i=0}^{1}u_{i}}{k})$ and $-d$, which are both independent of $i$, drop out of the quotients and we are left with showing that

[TABLE]

which is equivalent to

[TABLE]

Since $\mu$ is assumed to be increasing in $i$ this inequality holds. ∎

We will assume that the rate of change of the feedback variable $\varphi$ can be written as

[TABLE]

where the constant $\epsilon_{e}$ will later be assumed to be sufficiently small and $h$ is continuously differentiable. We make the following “negative feedback” assumption:

1. (NF1)

The function $h$ satisfies $h<0$ when $\varphi>M$, and $h>0$ when $\varphi<M$.

An example of such a function is $h(u)=M(u)-\varphi$, but it is not necessarily the case that $h$ can be expressed as a function of $M$. The NF1 assumption of course implies that for each value of $M$ there is a unique value $\varphi=\varphi(M)$ for which $h=0$. The assumption that this unique value equals $M$ is merely a convenience, which by the Implicit Function Theorem can be obtained by a change of coordinates whenever $\frac{\partial h}{\partial\varphi}(u,\varphi)\neq 0$ is satisfied for all $u$ and $\varphi=\varphi(M)$.

We will make the following additional assumptions guaranteeing both a negative and a positive feedback.

1. (NF2)

There exist $1>\varphi_{+}>\varphi_{-}>0$ with the following property: For any given non-zero $u$, the rates $\mu$ are decreasing in $i$ for $\varphi>\varphi_{+}$, and increasing in $i$ for $\varphi<\varphi_{-}$. 2. (PF)

For each $\varphi$ there is a unique value $M_{\varphi}$ such that $\mu$ is strictly increasing at $(u,\varphi)$ whenever $M(u)>M_{\varphi}$, and strictly decreasing whenever $M(u)<M_{\varphi}$. We assume that $M_{\varphi}$ is non-decreasing with $\varphi$.

In other words, when $\mu$ is non-decreasing, it must remain so when $\varphi$ is decreased or when $M(u)$ is increased, and in the latter case must become strictly increasing. Similarly, when $\mu$ is non-increasing, it must remain so when $\varphi$ is increased or when $M(u)$ is decreased, in the latter case it must again become strictly decreasing.

By continuity of $\mu$ it follows from (PF) that $\mu$ is constant when $M(u)=M_{\varphi}$. Note that we may redefine $\varphi_{+}$ as the smallest $\varphi$ for which $M_{\varphi}=1$, and similarly $\varphi_{-}$ as the largest $\varphi$ for which $M_{\varphi}=0$. Note also that we do not assume that $\mu$ is a function of $M$.

Let us now add mutations to the model:

[TABLE]

We assume that each function $g_{i}$ is non-negative when $u_{i}=0$, is decreasing in $u_{i}$, and is strictly increasing in each $u_{j}$ for $j\neq i$. In particular $g_{i}$ is strictly positive when $u_{i}=0$ but $u\neq 0$, thus mutations from any strain to any other strain are possible. We will later discuss alternative assumptions, making it possible to restrict some of the mutations.

Lemma 7**.**

Each $u_{i}$ remains bounded from below by a constant of order $O(\epsilon_{m})$.

Proof.

It follows from the formula for $\frac{du_{i}}{dt}$ that $\|u\|:=\sum_{i=0}^{1}u_{i}$ remains bounded from above and below, i.e. $\|u\|=O(1)$. It follows that as $u_{i}\rightarrow 0$:

[TABLE]

Thus for $u_{i}$ small, the worst case scenario is that

[TABLE]

for some uniform constants $c,C>0$. It follows that

[TABLE]

∎

Theorem 8**.**

For $\epsilon_{e}$ and $\epsilon_{m}$ sufficiently small there exist orbits $(u(t),\varphi(t)$) for which $M(t)$ fluctuates arbitrarily closely between [math] and $1$.

Proof.

Suppose that we start with initial values $u(t_{0}),\varphi(t_{0})$ for which $M(t_{0}):=M(u(t_{0}))>\varphi_{+}$, and for which $\mu$ is increasing in $i$. It follows that $\varphi(t_{0})<\varphi_{+}$, hence $\varphi(t_{0})<M(t_{0})$ and therefore $\varphi$ is increasing at time $t_{0}$.

By choosing $\epsilon_{m}$ sufficiently small we can guarantee that $M(t)$ remains arbitrarily close to $1$ until $\varphi$ is arbitrarily close to $\varphi_{+}$, say $\varphi>\varphi_{+}-\delta$ for $\delta>0$ arbitrarily small. Note that $\varphi$ remains increasing while this is the situation, hence at some time $t_{1}>t_{0}$ we must have $M(t_{1})=M_{\varphi(t_{1})}$. We may assume that $t_{1}$ is the first time that this occurs, from which it follows that $M(t_{1})$ must still be arbitrarily close to $1$, while $\varphi(t_{1})$ must be arbitrarily close to $\varphi_{+}$. In particular $\varphi$ is still increasing at time $t_{1}$, while $M(t)$ is decreasing due to mutations. As a consequence (PF) implies that $\mu$ becomes strictly decreasing, hence $M(t)$ will continue to decrease.

Since $\varphi$ remains increasing and $M(t)$ remains decreasing, it follows that there is a smallest time $t_{2}>t_{1}$ for which $M(t_{2})=\varphi(t_{2})$. We may assume that $t_{2}$ is the first time at which equality occurs, from which it follows that $\varphi$ has only increased between $t_{1}$ and $t_{2}$, and hence $\varphi(t_{2})>\varphi_{+}-\delta$.

We claim that $u_{0}(t_{2})$ is bounded from below by a constant that is independent from $\epsilon_{m}$. To see this, note that by the assumption that the vector $u$ remains bounded, it follows that the $u_{i}$’s for $i\neq 0$ must remain arbitrarily small for $t\in(t_{0},t_{1})$. It follows that the corresponding growth factors $\mu_{i}\left(1-\frac{\sum_{i=0}^{1}u_{i}}{k}\right)-d$ must remain strictly negative, with a uniform bound from above. It follows that populations $u_{i}$ for $i\neq 1$ remain comparable to $\epsilon_{m}$, where Lemma 7 implies the estimate from below, and in particular the populations $u_{i}(t_{1})$ for $i\neq 1$ are comparable to each other, with ratios independent of $\epsilon_{m}$. Recall that for $t\in[t_{1},t_{2}]$ we noted that $\mu$ is increasing, and hence the growth factor $\mu_{i}\left(1-\frac{\sum_{i=0}^{1}u_{i}}{k}\right)-d$ is largest for $i=0$. It follows that in the interval $[t_{1},t_{2}]$ the population $u_{0}$ grows faster than any other population $u_{i}$, with a strictly larger exponential coefficient. At time $t_{2}$ the average $M(t)$ has decreased by an amount independent of $\epsilon_{m}$. Since the total population remains bounded away from [math] by assumption, it follows that the size of $u_{0}$ must have increased by an amount independent of $\epsilon_{2}$, thus obtaining the claim.

By assuming that $\epsilon_{m}$ and $\epsilon_{e}$ are sufficiently small, it follows from (PF) and continuity of $\mu$ that $\mu$ will remain decreasing and $M(t)$ decreases below $\mu_{-}$, say at time $t_{3}$, and that the time interval $t_{3}-t_{2}$ is bounded and independent of $\epsilon_{m}$ or $\epsilon_{e}$. Since $\epsilon_{e}$ is assumed to be small, it follows that $\varphi(t_{2})\sim\varphi(t_{1})$. We have ended up with assumptions on $M(t_{3})$ and $\varphi(t_{3})$ that are symmetrical to those on $M(t_{0})$ and $\varphi(t_{0})$. The symmetry of our assumptions implies that the process will repeat itself, causing arbitrarily large fluctuations in $M(u)$. ∎

If we drop the assumption that mutations from any strain to any other strain are possible, and replace it instead by the much weaker assumption that given any two strains $u_{i}$ and $u_{j}$ there is a possible sequence of mutations from $u_{i}$ to $u_{j}$, Lemma 7 fails, and the above proof breaks down. We cannot guarantee that at time $t_{2}$ the strain $u_{0}$ has increased to a definite size, independent of $\epsilon_{m}$, and as a result we cannot give a bound on the time interval $t_{3}-t_{2}$.

This issue can be solved by assuming that the constant $\epsilon_{e}$ is sufficiently small, where the bound on $\epsilon_{e}$ may have to depend on the choice of $\epsilon_{m}$. Note the difference with the above statement, which holds whenever both $\epsilon_{m}$ and $\epsilon_{e}$ are sufficiently small. In practice the stronger assumption on $\epsilon_{e}$, which can imply that $\epsilon_{e}$ is much smaller than $\epsilon_{m}$, may or may not be desirable.

5. Examples: Rate equations, diversity index and supplementary figures

5.1. Shannon index

The Shannon diversity index (H) is used in ecological research to describe species richness. Here we used the measure to describe phenotype richness within a species (we do not discuss whether bimodal phenotype distributions lead to different species). We calculate this diversity with the following formula:

[TABLE]

Where $n$ is the number of phenotypes and $p_{i}$ the proportion of phenotype $i$ in the population.

5.2. One species with two phenotypes and a fixed population size

In section 4.1 we give a proof for a system with two phenotypes and a fixed population size, and in this section we made an example of such a system. For this example we used the following equations:

[TABLE]

$R$ and $\varphi$ are short for $R{t}$ and $\varphi{t}$ since they are time dependent. $R$ is the ratio of the two phenotypes, $A$ and $B$ ($r=\frac{A}{B})$, where $A$ is the phenotype interacting with the negative feedback $\varphi$. Figure 10 shows the result of time simulations of the example.

5.3. One species with two phenotypes and a variable population size

In section 4.2 we give a proof for a system with a variable population size. An example is shown in Figure 2 in the main text and the equations for the time derivatives of the populations with phenotypes $A$ and $B$ and the inhibitor $\varphi$ for that figure are:

[TABLE]

$A$, $B$ and $\varphi$ are all time dependent and therefore short for $A(t)$, $B(t)$ and $\varphi(t)$. $k_{A}$ and $k_{B}$ are parameters describing the growth of $A$, $d$ is a death rate and $K$ the carrying capacity. The parameters used for Figure 2 are $k_{A}=0.01$, $k_{B}=0.5$, $d=0.3$ and $K=10$. For the slow feedback and evolution we used $\epsilon_{E}=0.0005$ and $\epsilon_{M}=0.00005$, for the fast feedback and evolution we used $\epsilon_{E}=0.1$ and $\epsilon_{M}=0.01$ and for no evolution we used $\epsilon_{E}=0.01$ and $\epsilon_{M}=0$.

5.4. Multi-varied trait

The equations used for Figure 4 in the main text are:

[TABLE]

The collection of phenotypes $I$ consists of $n$ phenotypes with value $i$ and density $u_{i}$. The total population size is $U=\sum_{i}u_{i}$. The average value of the trait in the population is $M=\frac{\sum_{i\in I}i\cdot u_{i}}{U}$. $\varphi$ is the level of the feedback compound, $K$ the carrying capacity and $d$ the death rate.

The rate with which mutations change the phenotype decreases exponentially with the difference between the phenotypes. This is to simulate that mutations with small effect are more prevalent than mutations with large effects.

We did a time simulation with different values of the rate of evolution ($\epsilon_{M}$) and the delay in the slow feedback ($\epsilon_{E}$). For Figure 4 in the main text we used the following parameters: $K=1$ and $d=0.01$ throughout and for the the RQ dynamics $\epsilon_{E}=10^{-4}$ and $\epsilon_{M}=10^{-6}$, for the fast feedback and evolution $\epsilon_{E}=0.1$ and $\epsilon_{M}=0.01$ and when there are no mutations $\epsilon_{E}=10^{-4}$ and $\epsilon_{M}=0$. We assumed a population to go extinct if the population size was under 0.005.

5.5. Two feedbacks

We simulated a system with 2 traits associated with 2 feedbacks for Figure 5 in the main text. The time derivatives of the phenotypes $u_{i}$, where $i=\{i_{1},i_{A}\}$ is a vector of the two phenotypes, are:

[TABLE]

The collection of phenotypes $I$ consists of $n$ phenotypes with value $\{i_{1},i_{A}\}$ and density $u_{i}$. The total population size is $U=\sum_{i}u_{i}$. The average value of the first trait in the population is $M_{1}=\frac{\sum_{i\in I}i_{1}\cdot u_{i}}{U}$ and of the second trait $M_{A}=\frac{\sum_{i\in I}i_{A}\cdot u_{i}}{U}$. $\varphi_{1}$ is the level of the feedback compound for the first trait, $\varphi_{A}$ the level of the feedback compound for the second trait, $K$ the carrying capacity and $d$ the death rate.

Since mutations are rare we ignore mutations in both traits at the same time and the rate of mutations in one trait from one phenotype to the other decreases exponentially with increasing difference between the phenotype.

We used the parameters $d=0.01$, $K=1$, $\epsilon_{E}=0.0001$ and $\epsilon_{M}=0.0005$ for Figure 5 in the main text.

5.6. Literature models

We used the competitor-competitor model from [22] where they show RQ dynamics (Fig. 2 in [22]). We changed the model to allow for a polymorphic population, which leads to the following equations:

[TABLE]

The parameters are given in [22] Eqs 6 and Fig. 2. $\epsilon_{E}$ is set to 1 in Main Text Fig. 6B and to 0.1 in Fig. 6C. The trait values ($i$ and $j$) range from 0.5 to 1.5 and we simulated 100 different phenotypes per species within this range.