How long do Red Queen dynamics survive under genetic drift? A comparative analysis of evolutionary and eco-evolutionary models.

BACKGROUND: Red Queen dynamics are defined as long term co-evolutionary dynamics, often with oscillations of genotype abundances driven by fluctuating selection in host-parasite systems. Much of our current understanding of these dynamics is based on theoretical concepts explored in mathematical models that are mostly (i) deterministic, inferring an infinite population size and (ii) evolutionary, thus ecological interactions that change population sizes are excluded. Here, we recall the different mathematical approaches used in the current literature on Red Queen dynamics. We then compare models from game theory (evo) and classical theoretical ecology models (eco-evo), that are all derived from individual interactions and are thus intrinsically stochastic. We assess the influence of this stochasticity through the time to the first loss of a genotype within a host or parasite population.
RESULTS: The time until the first genotype is lost ("extinction time"), is shorter when ecological dynamics, in the form of a changing population size, is considered. Furthermore, when individuals compete only locally with other individuals extinction is even faster. On the other hand, evolutionary models with a fixed population size and competition on the scale of the whole population prolong extinction and therefore stabilise the oscillations. The stabilising properties of intra-specific competitions become stronger when population size is increased and the deterministic part of the dynamics gain influence. In general, the loss of genotype diversity can be counteracted with mutations (or recombination), which then allow the populations to recurrently undergo negative frequency-dependent selection dynamics and selective sweeps.
CONCLUSION: Although the models we investigated are equal in their biological motivation and interpretation, they have diverging mathematical properties both in the derived deterministic dynamics and the derived stochastic dynamics. We find that models that do not consider intraspecific competition and that include ecological dynamics by letting the population size vary, lose genotypes - and thus Red Queen oscillations - faster than models with competition and a fixed population size.

Background

Diversity, induced by continuous co-evolution can theoretically be maintained by the intense antagonistic relationship of hosts and parasites. This is the central part of the Red Queen hypothesis, verbally first formulated by van Valen in 1973 [1]. The hypothesis has been mathematically formulated in many models. However, owing to the modern usage of the term ‘Red Queen’ for different but related phenomena [2-9], the models have diverging foci and many lack the implementation of stochastic forces and ecological dynamics. A common synonym for the term Red Queen dynamics is fluctuating selection dynamics (FSD). Such fluctuations can be induced by co-evolving hosts and parasites and, as one possibility, be driven by negative frequency-dependent selection (NFDS), where host and parasite genotype abundances oscillate in time and every genotype can temporally be best adapted. In detail, since parasites are selected to target the most common resource, being a rare host genotype is advantageous. This temporary high fitness makes the genotype grow in relative abundance, but before it can take over the whole population, it is severely diminished by the profiting parasites genotypes, which target this now common host type. By contrast, in arms race dynamics (ARD) novel favoured genotypes spread in the entire population by recurrent selective sweeps. The terms NFDS and ARD are both referred to as Red Queen dynamics [10-12] and describe an ongoing co-evolutionary change without approaching an equilibrium. In this paper, we use the term Red Queen dynamics for NFDS, as is commonly done in the literature, but return to other definitions of the Red Queen in the discussion.

Although Red Queen dynamics is a well-known and frequently cited concept, there is only little evidence for the ubiquitous prevalence of long term Red Queen dynamics in nature – empirical challenges preclude the observation of more than a few subsequent oscillations, as these require a major amounts of intensive and challenging lab work [13-17]. Thus, most work on the actual long term temporal dynamics is theoretical, often dealing with evolutionary dynamics or epidemiological dynamics in a deterministic fashion. We have summarised some of the literature in the context of these assumptions in Table 1 (methods in the Additional file 1). Similar literature summaries exist with a focus on sexual vs. asexual reproduction [8] or host-parasite coevolution models [18]. Many theoretical studies build on evolutionary game theory [19] and a zero-sum assumption, where the harm done to the host equals the benefit for the parasite, which was already envisioned by van Valen at the time. Some of the models are implemented with equations that describe both species’ dynamics (explicit host-parasite HP dynamics), other studies, especially on the evolution and maintenance of sexual reproduction (Red Queen Hypothesis) revert to epidemiological models (susceptible-infected SI models), sometimes in the pursuit of including population dynamics. The present work focuses on evolutionary host-parasite models in comparison with eco-evolutionary models that include population dynamics without using SI models.

Table 1

Literature overview

Ref.	Authors (year)	focus	deterministic/ stochastic	equations/ method	population size
[31]	Schaffer and Rosenzweig (1978)	HP, CSS	deterministic	ODE	constrained4
[32]	Seger (1988)	HP, many genotypes, chaos	deterministic	RE	constant
[33]	Nee (1989)	HP, co-evolution, recombination	deterministic	RE	constant
[34]	Dybdahl and Lively (1998)	time lag, experiment	deterministic	RE	constant
[35]	Boots and Sasaki (1999)	infection on lattice	both	ODE, IBM, AD	variable
[36]	Peters and Lively (1999)	fluctuating epistasis	deterministic	RE	constant
[37]	Sasaki (2000)	multilocus GfG	deterministic	ODE	constant
[38]	Agrawal and Lively (2001)	HP, selfing vs outcrossing	deterministic	RE	constant
[39]	Agrawal and Lively (2002)	HP, GfG vs MA	deterministic	RE	constant
[40]	Gandon (2002)	HP, local adaptation (spatial)	deterministic	RE	constant
[41]	Gandon (2004)	SI, multihost parasites	deterministic	ODE, AD	variable
[20]	Kouyos et al. (2007)	HP, oscillations in stochastic model	both7	ODE	constant5
[42]	Alizon and van Baalen (2008)	multiple infections, within-host and SI	deterministic	ODE, AD	variable
[43]	Agrawal (2009)	HP, sex vs recombination	deterministic	RE	constant
[44]	Best et al. (2009)	SI, transmission, susceptibility	deterministic	ODE, AD	constant
[21]	Engelstädter and Bonhoeffer (2009)	HP, RQ oscillations	deterministic	RE	constant
[45]	Lively (2010)	sex (long term persistence)	both6	RE	variable
[46]	Greischar and Lively (2011)	HP, extinction risk	deterministic	RE	constrained
[47]	Gilman et al. (2012)	HP, multiple host traits, resistance	stochastic	IBM	constant, constrained4
[48]	Mostowy and Engelstädter (2012)	interaction matrices, sex, LD	deterministic	RE	constant
[28]	Gokhale et al. (2013)	HP, population size	stochastic	IBM	variable, constrained
[49]	Luijckx et al. (2013)	MA, Daphnia	deterministic	RE	constant
[50]	Abou Chakra et al. (2014)	HP, plastic behaviour	both	ODE, IBM	constant
[51]	Taylor et al. (2014)	HP, virus of virus	deterministic	ODE	constrained
[23]	Ashby and Gupta (2014)	SI, state-dependent sex, MA	deterministic	ODE	variable
[8]	Ashby and King (2015)	SI, diversity, transmission, sex	stochastic	IBM	constant
[52]	Engelstädter (2015)	HP, infection matrices	deterministic	RE	constant
[53]	Rabajante et al. (2015)	HP, many types	deterministic	ODE	constrained
[25]	Song et al. (2015)	HP, population size, GfG MA	deterministic	ODE	constant, variable
[54]	Hesse et al. (2015)	environment, specialisation	deterministic	ODE, AD	constrained4
[24]	Gómez et al. (2015)	oscillation vs. arms race	stochastic	IBM	variable
[55]	Rabajante et al. (2016)	HP, rare types	deterministic, noise1	ODE, SDE	constrained
[56]	Nordbotten and Stenseth (2016)	HP, RQ vs stasis	deterministic	PDE	constrained4
[57]	Best et al. (2017)	SI, no specificity, FSD	deterministic3	ODE, AD	constrained4
[58]	Bonachela et al. (2017)	crossfeeding	deterministic2	ODE	constrained
[59]	Greenspoon and Mideo (2017)	relatedness, transmission	deterministic	RE	constant
[60]	Lively (2017)	allopatric, sympatric parasites	deterministic2	RE	constrained
[61]	Nuismer (2017)	local, global adaptation	deterministic2	RE	constant
[62]	Veller et al. (2017)	HP, speed of evolution (RQ, RK)	stochastic	IBM	constant
[63]	Ashby and Boots (2017)	HP, SI, GfG MA	deterministic	ODE	constrained4
[27]	MacPherson and Otto (2018)	SI, HP, MA, RQ oscillations	deterministic	ODE	constant, constrained4
[18]	Ashby et al. (2019)	HP, population size change	deterministic	ODE, AD	constrained
[]	Current paper	(HP, MA, RQ) population size, extinction time	stochastic	IBM	constant, constrained, variable

Mathematical models and properties discussed in this paper sorted by publication year. Many models deal with relative allele or genotype abundances without considering ecological dynamics – these have been categorised as constant population size models. Those models that do include a changing population size and stochastic effects mostly do not analyse the stability of long term oscillations which is the focus of this paper. (See the notes on this literature survey in the Additional file 1).

ODE/PDE/SDE: ordinary/partial/stochastic differential equation, IBM: individual based model (stochastic simulations), RE: recursion equation, SI: susceptible-infected (epidemiological) model, HP: explicit host-parasite model, AD: adaptive dynamics (most often ODE with added mutants), MA: matching alleles, GfG: gene for gene, RQ: Red Queen (oscillations in genotype abundances or in trait space), RK: Red King (slow evolution favoured), CSS: coevolutionary stable strategy.

1not intrinsic stochasticity

2stochastic mutants added

3adaptive dynamics simulations (no intrinsic stochasticity)

4via carrying capacity (density dependent death or competition term)

5but discussed

6some randomness in infection (±1 in next generation)

7when time discrete, only host stochastic

While many studies assess the occurrence of oscillating selection dynamics and show under what assumptions oscillations dominate [18, 20–27], only few studies include both ecological population dynamics and stochastic noise, although the combination of the two has been shown to result in a fast loss of genotypes in either population [28]. It has been difficult to derive a stochastic model that easily switches between constant and changing population size using a single parameter. For example Gokhale et al. [28] artificially normalised population size every few generations. Here, we take a different approach and compare the modelling framework of evolutionary game theory, where population size is constant by design, to eco-evolutionary dynamics from the field of theoretical ecology, where population size is inherently free to change over time. Our goal is not to present the one model that is the best description of reality, but to illustrate how different modelling assumptions can drive the results from such models.

Specifically, we use individual-based models, since ecological and evolutionary dynamics of populations are driven by events on the individual level. The models are based on haploid and asexual populations that live in a well-mixed environment where encounters are density dependent. Individuals are born, interact with other individuals of their own or opposing species and die. Generally, we will consider at least two genotypes and track the associated abundances H1,H2,P1,P2 and the total population sizes N,N of hosts and parasites over time. We simulate the dynamics using a uniformly distributed initial standing genetic variation and the simple matching allele interaction profile, where parasites are highly specialised [29, 30] on a particular host genotype and identical in all other aspects. Yet, the way this interaction profile enters in the dynamical equations and thus defines fitness for the individual genotypes is very different between the models. In population dynamics models these events happen at constant rates and depending on the density of the interacting individuals. A similarly simple, yet completely different approach is the stochastic birth-death process which tracks only the evolutionary dynamics. In each time step one individual is born, proportional to its current ‘fitness’ and another individual dies proportional to its density.

These models all produce Red Queen dynamics (NFDS) and we assess the robustness of those by measuring the time to extinction, which we define as the earliest time that any genotype from the initial genetic variation is lost in either population. Further, we consider the impact of the derived deterministic dynamics and the influence of ecology in the form of a population-size-change on this extinction time. The time to extinction of a genotype represents the durability of the stochastic oscillations. Without the immigration or re-emergence of extinct genotypes, the diversity of both populations declines in the long run.

Results

Evolutionary dynamics depict the change of relative genotype abundances over time and can be examined without keeping track of population size changes. However, it is well known that ecological dynamics can feed back on evolutionary dynamics. We want to understand this feedback in the context of Red Queen dynamics. To this end we compare models from evolutionary game theory, that do not include population size changes and theoretical ecology models that do. The models have been widely used in the literature and represent the simplest case of Red Queen dynamics with a matching allele interaction profile (for details see the “Methods” section below and Additional file 1).

The matching-allele host-parasite Red Queen dynamics in evolutionary and eco-evolutionary models

In an evolutionary birth-death process one individual is born and another dies in each population, here host or parasite, and in each time step. Thereby, population size remains constant and the focus lies on the genotypic composition of a population. The Evo + and Evo processes (see Table 2 and in the “Methods” section for a definition) are such birth-death processes [64, and references therein; 65]. Individuals are chosen to die randomly, but the individual that reproduces is chosen proportional to the fitness advantages of that genotype relative to other genotypes in the population. The fitness effects are imposed by the current state of the antagonist population and an interaction matrix. In the Evo + process, the fitness effect is normalised by the average fitness effect over the whole host population, which leads to a kind of intra-specific competition (+) while in the Evo process the difference in fitness effects is compared between a genotype-one individual and a genotype-two individual, thus competition is pairwise. Because of the population size constraint, both processes can be analytically treated (see Additional file 1) when implemented in discrete time (prefix dt).

Table 2

Model overview

Model	description	features	stochastic dynamics (small N) ×	deterministic dynamics (large N) ∗	stochastic dynamics (medium N) †	population size change
evolution. game theory
Evo +	Birth-death process. Which individual reproduces depends on the current fitness effect by the antagonist, normalised by the population average fitness effects (+)	intraspecific competition (+)	slow extinction	stasis	NFDS	no
Evo	Like Evo + but fitness effects are compared between two individuals not with the population average	pairwise competition	slow extinction	NFDS	extinction	no
Hybrid	Hybrid model with reactions between two genotypes of different populations, single birth of parasite and death of host by dynamically adjusted rates.	no competition	extinction	NFDS	extinction	yes, but dynamically constrained
theoret. ecology
EcoEvo +	Independent reactions between individual hosts and parasites, single birth and death events or competition in hosts	intra-host competition (+)	fast extinction	stasis	NFDS	yes, but carrying capacity
EcoEvo	Like EcoEvo + but without competition within hosts. For infinite population size this is the Lotka-Volterra dynamics	no competition	fast extinction	NFDS	extinction	yes, uncon-strained

Model names and their main differences. The Evo + and Evo model are derived from evolutionary game theory while the EcoEvo + and EcoEvo model stem from theoretical ecology. The Hybrid model combines elements from both. Models are ordered by population size constraint. The deterministic dynamics apply to the two types matching alleles interaction matrix. Details on the models and analysis are available in the Additional file 1.

×population size change speeds up the extinction of genotypes (Fig. 1)

∗for large population sizes N the deterministic dynamics dominate (Fig. 2). Damped oscillations lead to an attractive equilibrium (‘stasis’). When the equilibirum is neutral genotype abundances oscillate induced by negative frequency-dependent selection (‘NFDS’).

†when population size is intermediate dynamics are strongly influenced by their deterministic characteristics but with stochastic noise. Stochastic dynamics with oscillations are stabilised by the attractive deterministic fixed point which can countervail the stochastic outward pull, postponing extinction (‘NFDS’). Without the attractive pull the time to the first extinction of a genotype is much shorter (‘extinction’).

Fig. 1

Example run illustrating that extinction is faster with ecological dynamics. Oscillations of host and parasite genotype abundances in the Evo + process with constant population size and EcoEvo + process with changing population size. The simulations start with an equal abundance of both genotypes H1(0)=H2(0)=N/2 and P1(0)=P2(0)=N/2. Method: Simulation of the stochastic processes with the Gillespie algorithm. Parameters: Total population sizes N=50,N=150 (only initially for the EcoEvo + model), selection strengths w=0.5,w=1, matching allele parameters α=1,β=0, death rate of the parasite d=1, birth rate of the host b=6, carrying capacity K=100, interaction rate , intra-specific competition rate . See the “Methods” section and Additional file 1 for method and parameter details

Fig. 2

Large population size limit. Relative abundances of two genotypes of host h1 and h2 and parasite p1 and p2 over time (left) and 2D representation (right) in the deterministic equivalents of the Evo + and Evo process with constant population size. Top: Intraspecific competition within the whole population (+) results in an attracting fixed point which is reached eventually and does not changed once reached, leading to stasis (also EcoEvo +). Bottom: Pairwise competition between individuals allows for a neutrally stable fixed point which neither attracts nor repulses the dynamics resulting in continuous co-evolution in the form of negative frequency-dependent selection dynamics (NFDS) around the internal fixed point (also EcoEvo). Method: integration of ordinary differential equations, the adjusted replicator dynamics (Evo +) and the replicator dynamics (Evo), which are the deterministic limits of the respective stochastic processes. Parameters: selection strength w=w=1, matching allele parameters α=1 and β=0

In models adapted from theoretical ecology the events of birth, death and interaction happen independently with external rates and, importantly, between populations (EcoEvo, comparable with the Lotka-Volterra dynamics in [28]). Host and parasite individuals encounter one another based on their densities and if they match, an interaction is carried out with a constant rate upon which a host dies or a parasite reproduces. When competition between hosts (+) is included, the host population grows logistically with a carrying capacity K. The host population size N reaches the carrying capacity in the absence of parasites. However, in the presence of parasites the population size is smaller than K because of the additional mortality from parasites (see Table 3 for details on the parameters in all models).

Table 3

Model parameters

Parameter	interpretation	default value / range	Models
Hi	absolute abundances of host genotype i	∈[0,NH]	all
Pi	absolute abundances of parasite genotype i	∈[0,NP]	all
hi	relative abundances of host genotype i	∈[0,1]	all
pi	relative abundances of parasite genotype i	∈[0,1]	all
NH	total host population size ΣiHi	-	all
NP	total parasite population size ΣiPi	-	all
bH	host birth rate	-	EcoEvo +, EcoEvo
dP	parasite death rate	1	EcoEvo +, EcoEvo
λ	rate of host death or parasite birth when genotypes match	-	EcoEvo +, EcoEvo
μ	intraspecific competition rate	-	EcoEvo +
K	carrying capacity (host population size in the absence of parasites)	-	EcoEvo +
α	fitness gain (loss) of matching genotypes for parasite (host)	1	Evo +, Evo, Hybrid
β	fitness gain (loss) of mismatching genotypes for parasite (host)	0	Evo +, Evo, Hybrid
wH,wP	selection intensity	∈[0,1]	Evo +, Evo, Hybrid
π_H_i,π_P_i	payoff from the game, defined for each genotype	-	Evo +, Evo, Hybrid
f_H_i,f_P_i	‘fitness’ from the game, defined for each genotype	-	Evo +, Evo, Hybrid
f_H¯,f_P¯	average ‘fitness’, defined for each population	-	Evo +, Evo, Hybrid
d_H_i	death rate of host genotype i	-	Hybrid
b_P_i	birth rate of parasite genotype i	-	Hybrid

Both Evo and EcoEvo modelling approaches are combined in an intermediate model with self-controlled, but not fixed, population size (Hybrid). The model is implemented as an individual interactions model, where reactions take place also between populations, but the rates of these events are taken from the game theory models: Host death and parasite birth happen according to the fitness effects, host death and parasite birth rates are then adapted dynamically to keep population size nearly constant.

From the derivations of the models (details in the “Methods” section and Additional file 1) some basic properties of the dynamics are obtained and summarised in Table 2. The evolutionary game theory models have a constant population size by design, whereas population size can change in all other models. The average behaviour of the individual-based stochastic processes is captured in the deterministic selection term and the noise term, which together determine the stochastic dynamics. The noise term is discussed in the “Results” section (Fig. 1). The role of intra-specific competition in the deterministic part is discussed in the “Discussion” section (Fig. 2).

In models with ecological dynamics genotypes die out faster

It is clear that populations with low total population sizes are more prone to genetic drift and the loss of genotypes than large populations. We now show that it is not only the population size but the possibility of population size change that speeds up this process. We exemplify our argument here, but provide a more detailed analysis below. As an example with two host genotypes and two parasite genotypes we select the Evo + process and the EcoEvo + process (Fig. 1). Starting with an equal abundance of genotypes we measure the time to the first loss of a genotype. When a genotype has died the population becomes monoclonal and oscillations are no longer possible. With the fixed population size in the Evo + process oscillations survive longer than in the EcoEvo + process with a changing population size. The evolutionary dynamics are similar and defined through the relative abundance of the types, but the population size change can speed up the frequency of event occurrences and increase the probability of extinction through the bottleneck effect when population sizes are low.

Intraspecific competition stabilises negative frequency-dependent selection

The equations that define the stochastic process consist of a deterministic selection term and a noise term and represent the mean and variance of many individual simulations. Therefore, it is impossible to understand the stochastic model without making the deterministic dynamics clear. Furthermore, when population size is large, the stochastic process approaches the more manageable deterministic dynamics (details in the Additional file 1). The deterministic equations for all models from Table 2 have an internal co-existence fixed point, where both genotypes exist in a fixed ratio, which does not change over time. This point is only attractive, if starting with suitable initial compositions of genotypes the dynamics approach the state, in this case in the form of damped oscillations. The intraspecific competition (+) in the Evo + process and the EcoEvo + process result in such an attractive pull (Fig. 2). A second possibility is neutral stability, where genotype abundances oscillate with a constant amplitude and period, which depend on the initial abundances. These neutral cycles are produced by models where individuals only compete with other genotypes locally like in the Evo process with pairwise competition or the EcoEvo model with no intraspecific competition and the Hybrid model (Fig. 2).

In our models, the noise in the stochastic dynamics always leads to extinction, while deterministic dynamics never do. When population size is large enough to be impacted by the deterministic behaviour but stochastic noise still plays a role, the global competition models (+) show persisting Red Queen oscillations. The deterministic ‘pull’ and the stochastic ‘push’ balance [66], prolonging extinction times. For models with neutral oscillations (NFDS) in the deterministic dynamics stochastic effects will on average increase the amplitudes and push the trajectories to the edges of the space towards a faster extinction of genotypes.

The single simulations (Fig. 1) are only a snapshot and one specific realisation of the stochastic processes. Ideally, we would analytically derive general extinction times depending on the parameters of the model. Yet, to derive an exact analytical solution for this problem is extremely challenging. In addition to simulations, we have calculated the numerical (but exact) extinction times for low population sizes and provide an approximative method based on the averaged noise (see Additional file 1 for further details). These methods are limited to a subset of the models and can thus not be used for a comparison of all models, but only to support the computationally costly simulations which provide the now following main result.

The strength of random effects depends on the model properties

We simulate 1000 replicates for a set of parameter combinations of the models with two genotypes in each population and record the time it takes until one genotype has died out. As a general trend, the more constrained a population size is, the longer oscillations survive (higher extinction times in Fig. 3). This holds true for small to intermediate population sizes – note a similar vertical order of extinction times to the ordering of models by population size constraint in Table 2. When population sizes become larger and the deterministic model properties gain influence, models with competition terms (+, stasis, compare with Fig. 2) have higher extinction times and therefore more stable Red Queen oscillations.

Fig. 3

Extinction time of either genotype of either host or parasite population for different initial population sizes of the parasite N for all models. We show the mean extinction time of any genotype over 1000 independent simulations (fat dots) and the distribution of those extinction times (shaded histogram area around the mean). The simulations start with equal abundance of both genotypes H1(0)=H2(0)=N/2 and P1(0)=P2(0)=N/2. Lines denote approximate results based on the average noise (see Additional file 1). The discrete time processes are simulated for values of N for which analytical results are valid. The Evo + process is not simulated for high parasite population sizes since the computation becomes extremely time-consuming and the trend is already clear. Parameters as in Fig. 1 except N=250,K=500, birthrate b∈{0.24,0.32,...,1.6} in the EcoEvo + and for the EcoEvo model b∈{0.12,0.16,...,0.8} and μ=0 to achieve the population sizes N displayed

By design, the discrete time (dt) processes have much higher extinction times and are thus not directly comparable to the continuous time simulations. A scaling would be possible for equal population sizes, but with different extinction routes and N≠N no such factor can be derived. The dtEvo + and dtEvo extinction times in Fig. 3 can therefore only be compared between them. For growing N, the dtEvo + process has an increased extinction time because of the stabilising attractive fixed point. This trend is even more pronounced in the approximate analytic solution (solid lines), inspired by Claussen [67, 68] (see Additional file 1). The error of the analytical approach cannot be neglected, but the qualitative trend is clearly visible and the result is fully analytical.

Due to the challenges of employing an exact analytical approach, we cannot perfectly tune the models for the same amplitudes, fluctuations and frequencies/periods of oscillations. The specific choice of the parameters is not necessarily directly comparable, but we have made an effort to choose them in a meaningful way, such that the fixed points are exactly the same and amplitudes comparable. We choose strong selection for the parasite w=1 and weaker selection for the host w=0.5 in the models derived from game theory, because the EcoEvo + model is built in a similar way: Parasite birth can only occur through the antagonistic interaction, but host mortality is also influenced by the competition term. While the parasite is obligate and thus completely dependent on the host, the host suffers, but does not always die from an infection.

The impact of selection intensities on the extinction times is further explored in the Additional file 1. We find that strongly diverging host and parasite selection intensities can counter-intuitively lead to more stable dynamics in the Evo process than in the Evo + process.

Diversity inflow results in sequential negative frequency-dependent selection dynamics and arms race dynamics

So far we have compared models with two genotypes in each species. We now provide an outlook of how diversity changes for many genotypes. We simulate one possible example with an initial uniform distribution of twenty genotypes in each species (see Additional file 1). Diversity, simply defined as the number of genotypes present in the population, declines exponentially with time at a constant rate. The manual re-introduction of an extinct, but temporarily best adapted parasite genotype can result in a selective sweep that leaves the parasite population monoclonal.

In reality, our genotypes are not as static in their traits as described here, but one of our ‘genotypes’ can actually be seen as an average of several individuals with slightly different traits. We now add a form of mutation or recombination to the model so that reproduction does not necessarily result in a clonal daughter, but a new individual with different traits. For example, parasites could evolve quickly by allowing beneficial mutations to produce other, even extinct, genotypes. Depending on the model system, a sexually reproducing host could also store genetic material to revive long extinct phenotypes by recombination. We abstract this by inserting a conversion rate μ from one genotype to the neighbouring genotype. For example with five pre-defined genotypes we have and and so on. The dynamics we observe now (Fig. 4) are not pure negative frequency-dependent selection dynamics, but a mixture of oscillations and arms race dynamics, where selective sweeps can make a population monoclonal in a very short time, but a re-introduction of extinct genotypes allows for oscillations to re-emerge.

Fig. 4

Negative frequency-dependent selection and arms race dynamics. Revival of genotypes and evolution of host (top) and parasite (middle) populations with five possible genotypes each. With the rate μ=0.005 and μ=0.01 genotypes convert to neighbouring genotypes through mutation or recombination. Stacked plots – evolutionary dynamics: the area covered by one colour is proportional to the relative abundance of that genotype of host or parasite at that time. Lower panel – ecological dynamics: total abundance of hosts and parasites. Method: the example is a stochastic simulation (Gillespie algorithm) of an EcoEvo + process. The simulations start with equal abundance of all five genotypes H(0)=N/5 and P(0)=N/5 for i=1,2,3,4,5. Parameters: N=300,N=900 (both initially), b=6,d=1,K=600,λ0=10

Discussion

In this paper we compare evolutionary models from evolutionary game theory to eco-evolutionary models from theoretical ecology to understand the impact of ecology and other model properties on the long term co-evolutionary Red Queen oscillations of host and parasite genotypes. The models are individual-based and intrinsically stochastic, thereby allowing genetic drift and the loss of genotypes from a population. Starting with an initially uniform distribution of genotypes, we define the extinction time as the first time that any genotype is lost from any of the two populations, and use this extinction time to measure the robustness of the Red Queen cycles and therefore, the maintenance of diversity. Our main result is that including ecology in models, in the form of a changing population size, leads to a faster loss of genotypes, when stochastic dynamics are considered. This result is similar to the simulation results by Gokhale et al. [28], where ecological dynamics were artificially removed from the simulations, in an attempt to make a straightforward comparison of eco-evo and evo dynamics. In contrast, we compare two modelling frameworks with historically developed differences between them. The models presented here are all based on the same widely used biological assumptions – haploid well-mixed host and parasite genotypes that interact through the matching-alleles infection matrix – but with differences in their mathematical properties: discrete and continuous time models with attractive or neutral deterministic dynamics. The models are further intrinsically stochastic, since they are derived from interactions between individuals. This inherent noise, genetic drift, also impacts the models within a given framework. The mean outward pull by noise that increases amplitudes and thus makes extinction more probable can be counteracted by intraspecific competition that pulls the dynamics back, decreases amplitudes and thus stabilises negative frequency-dependent selection dynamics, resulting in longer extinction times. Finally, we provide a snapshot of what happens when standing genetic variation is large initially. If no inflow of genotypes via mutation or migration is provided, the number of genotypes in a population will decline exponentially. However, when conversions between neighbouring types are allowed with a small mutation rate, negative frequency-dependent selection dynamics and selective sweeps can occur sequentially.

Previous theoretical studies have similarly examined the persistence of Red Queen oscillations under ecological feedbacks. For example Goméz et al. [24] found fluctuating selection and arms race dynamics in an epidemic model (host-focussed) with explicitly modelled parasite populations. MacPherson and Otto [27] also combined epidemiological and neutrally stable host-parasite dynamics and showed that this can dampen allele frequency oscillations, which leads to stasis in their deterministic model but would return to oscillations under stochasticity (see Table 2). Recently, the game theoretical replicator dynamics were mathematically tuned for population size influence using a single parameter [18] resulting in damped oscillations and thus stable polymorphism for both matching-alleles and gene-for-gene infection matrices. While we argue that eco-evolutionary feedbacks increase oscillation amplitudes, Ashby et al. argue that oscillation amplitudes are decreased over time. However, the population dynamics in their model were dampened by a maximal value which resembles our intra-specific competition leading to stasis (Fig. 2). Furthermore, their models are deterministic which closer resembles our models when population size is large, where stabilising effects have a larger influence. In the more theoretical literature, it is now well established that assumptions such as population size fluctuations and stochasticity can result in more rapid extinction ([69, 70], and many more).

The stabilising property of intra-specific competition is documented in the literature [71]. Intra-specific competition (+) enters in our evolutionary models as part of a genotype’s fitness effect that is compared to the focal population’s average fitness, whereas in the eco-evolutionary models it is implemented as an ecological intra-specific competition term. Both the evolutionary and the ecological implementation of this intra-specific competition stabilise the dynamics and lead to stasis following damped oscillations. The more commonly used host-parasite co-evolution models result in neutrally stable oscillations whereas damped oscillations are often seen as a termination of Red Queen dynamics. Yet, exactly this stasis shows similar oscillation patterns when stochasticity perturbs dynamics away from the stable fixed point (noise induced oscillations [66]). In a stochastic world, pure host-parasite dynamics therefore result in fast extinction, which would only be stabilised by intraspecific competition. For larger population sizes, when the stability of the fixed point gains in importance, the dynamics are pulled more towards the inner equilibrium state, making stochasticity less influential. Thus, only for organisms with large population sizes and good mixing, intraspecific competition would not be necessary for sustained Red Queen oscillations.

Although this study does not explicitly analyse modes of reproduction, our final result shows how reviving extinct genotypes can restore Red Queen dynamics. If parasites can evolve more quickly due to shorter generation times and larger numbers, then hosts are given an advantage by being able to “store” genotypes through recombination. Also, if clonal reproduction accumulates mutations (Muller’s ratchet), this could impact population sizes and sexual reproduction would be even more important [72, 73]. Ashby and King [8] devised a stochastic individual based susceptible-infected model with diploid sexual hosts and showed that high diversity cannot maintain sexual reproduction when parasite transmission rates are low. Although our models are more abstract concerning reproduction, we do explicitly model parasites. If parasite populations are well mixed and diverse with high mutation rates, this can again select for higher diversity through sex, like in [24], where fluctuating selection dynamics, and thus high diversity, is more likely when hosts encounter a diverse parasite population and the disease load is high. Furthermore, our models can include global competition in both species or resource competition in the host, which stabilise the oscillating dynamics. More support for recombination during parasite infection was shown in [48], where hosts could optionally switch between two modes of reproduction. See also [4] for a comprehensive connection to the Red Queen Hypothesis for sexual reproduction.

We have shown that in the same setting and with the exact same parameters sequential occurrences of oscillating selection and arms race dynamics are possible. We show only a snapshot and we do not quantify dynamics as is done in [24], but we find it to be an interesting aspect that the dynamics can occur temporarily in the same simulation, with the same settings and assumptions. The more complete picture could include all possibilities discussed in the Red Queen literature: there can be constant extinction, as suggested by van Valen on a taxonomic level and there can be oscillations and arms race dynamics as suggested by host-parasite interactions and the resulting co-evolution. With our preliminary results we might be going too far if we also justify sexual reproduction, yet, without recombination or mutation, diversity decline is inevitable theoretically.

Our models explore stochasticity under different restrictions of population size, while other modelling aspects are kept relatively plain. In the present work, the infection pattern is restricted to the matching alleles model, and the zero-sum assumption, yet this is necessary for oscillations [21]. Other infectivity patterns that result in Red Queen dynamics have not been examined here. The gene-for-gene infection matrix could show similar results since the oscillations are also neutrally stable, yet including ecological dynamics changes the complexity of the cycles [25]. In general the robustness of those cycles under stochasticity would depend on the details of how different infection mechanisms are translated into mathematical equations. Further limitations are the haploidy of both hosts and parasites and thereby asexual reproduction, the lack of life history or infection history and there is no spatial structure and evolution in the values of resistance or infectiousness. We do, however, briefly explore the effects of including more genotypes and mutation as a means to revive genotypes. There is an increasing effort to openly discuss how verbal models and biological assumptions enter into models [27, 74]. Making the assumptions clear and readily available should be the standard for future publications. For stochastic processes the analogous deterministic dynamics should be stated to provide the reader with a more complete picture of stochastic dynamics.

The model predictions presented here although quite abstract may nevertheless apply to the real world. Bottlenecks are likely more common in natural host-parasite associations [75] than usually assumed and, therefore, the interaction dynamics are likely shaped by genetic drift and, thus, stochastic effects. Eco-evolutionary feedbacks have been confirmed to impact the form of co-evolution in bacteria-virus experiments [76]. Increasing diversity in the parasite or higher exposure lead to a shift from negative frequency-dependent selection to arms race dynamics in two bacteria-phage systems [24, 77]. Oscillations alongside incomplete selective sweeps were recently even documented in a nematode-bacteria interaction [17]. It would now be of particular interest to assess the occurrence of bottlenecks, drift and competition in natural host-parasite associations and relate them to the resulting allele frequency dynamics. Such empirical data would help us to obtain a more general understanding of host-parasite co-evolution and potentially question the importance of sustained Red Queen oscillations in this context.

Conclusions

We have shown that models equal in their verbal biological description can be quite different in their mathematical details, with great consequences for both deterministic and stochastic dynamics. The loss of genotypes is inevitable in stochastic models without mutation or immigration. This extinction is faster when ecological dynamics are considered in an evolutionary model. Competition between genotypes within a species stabilises the dynamics and slows down extinction thus sustaining Red Queen dynamics. The applicability of models to the real world thus depends on the system of interest and mathematical details should be carefully considered for each particular case study. When bottlenecks, drift and competition are observed, the model needs to be adapted accordingly.

Methods

The following method descriptions are short explanations of the stochastic processes used in this manuscript. The precise equations and methods of analysis can be found in the Additional file 1. The simulation code is provided at https://github.com/HannaSchenk/ShortLifeRQ.

The discrete time Evo process (dtEvo), also discrete time Moran process, is a stochastic birth-death process, with a constant population size [78], often used in evolutionary game theory (see for example [79, 80] or [81]). Each birth-death reaction has a reaction probability (or transition probability), depending on the state of the system in each discrete time step Δt=1. The original definition ensures that the probabilities sum up to one so that one reaction (also reactions where no transition happens - when birth and death event happen within the same genotype) takes place in each time step. In the Moran process, a ‘payoff’ π is what a genotype gains from interactions with others. The interaction matrix is

where α is the fitness gain for the parasite and the fitness loss for a host if the genotypes match, and β is the fitness gain or loss for mismatching pairs (here α=1 and β=0). For example, the probability of a P1 birth and thus a subsequent P2 death is proportional to , where h1 and h2,p1 and p2 are relative abundances. How much this effects the so-called ‘fitness’ f is controlled with the selection intensity w (or w for the host) so that . This per capita ‘fitness’ is then normalised by a dynamically changing population average (this is the intraspecific competition step) and multiplied with the current abundance of the genotype. Thus is then the birth probability of a genotype one parasite. The death probability is simply density dependent, thus the total probability of replacing a genotype-two (death) by a genotype-one (birth) parasite is . Since we are modelling two populations (host and parasite), we choose to update both populations simultaneously instead of sequentially, such that frequency changes of host genotypes and parasite genotypes can happen at once. The deterministic limit (population sizes N,N→∞ and time steps Δt→0) of the Moran process in a single population is usually the differential equation of the replicator dynamics, however, in a two-population model the average fitness within each population is different and thus the adjusted replicator dynamics become the deterministic analogue [64, 82]. The adjusted replicator dynamics for host-parasite interactions have a globally attractive inner fixed point, in the symmetric matching alleles case this is the equal abundance of all genotypes.

The discrete time Evo process (dtEvo) [64, called pairwise comparison process or local update process in evolutionary game theory], is another birth death process, nearly equivalent to the Moran process but here competition is strictly local and pairwise, not normalised by a global average fitness. What is in the Moran process is here . The ‘fitness’ of parasite 1 only depends on the difference in ‘fitness’ to parasite 2 which depends on the abundances of host genotypes (see equation for ), but not, as when normalising with , on the relative abundances of the parasite genotypes. Thus, in the pairwise comparison process, the antagonist influences globally (since there is no spatial structure), but within a species the competition is local. This results in the recovery of the replicator dynamics with neutral cycles in the deterministic limit.

A Gillespie algorithm [83] can be employed to simulate the above stochastic processes. In this case reaction rates (not probabilities) are calculated for each species and, using random numbers, the shortest waiting time for each reaction is determined under the assumption of exponential waiting times. The reaction with the shortest time takes place and time is updated accordingly. This makes time continuous and time steps unequal. In contrast to the discrete time models, the Gillespie algorithm only updates one species at a time. The now following processes are also implemented using a Gillespie algorithm.

The EcoEvo process uses independent reactions of host birth, parasite death and host-parasite interactions similar to the individual-based equivalent of the Lotka Volterra equations. This results in a microscopic process often believed to be a more natural approach because the reactions describe individual and independent events on the ‘microscopic’ level rather than population dynamics on the ‘macroscopic’ level. Host birth reactions are density dependent with constant rate b, parasite death is density dependent with constant rate d and a density-dependent interaction of matching host-parasite pairs can result in the death of the host or the birth of a parasite with constant rate λ. The population size has no restrictions in this case and freely follows the evolutionary dynamics. The deterministic analogue has a neutrally stable fixed point.

The EcoEvo process are like the EcoEvo independent reactions, but with additional competition in the host population. Density-dependent interactions of two host individuals, independent of the genotype, result in the death of one of the individuals with constant rate μ. This model, when reduced to the deterministic limit, is an antagonistic interaction model with logistic growth in the host (carrying capacity K=b/μ) and an attractive inner fixed point.

The Hybrid model is a process with self-controlled population size. It is built on the EcoEvo model but with constrained birth and death rates adapted from the Evo model and dynamically varied to balance birth and death events on average. Thus, the population size is tightly controlled, yet it is not strictly constant. Building on the stochastic processes from evolutionary game theory above, one can set up a process that utilises the infection matrix for death events in the host and birth events in the parasite in explicit individual reactions. For example the dynamic reaction rate of a death event of a host genotype one is . The birth rates for the host and the death rates for the parasite are then dynamically adjusted to equal the total rates of host death and parasite birth. The deterministic limit is the replicator dynamics with a neutrally stable fixed point, as in the pairwise comparison process.

Additional file 1 Supplementary files. The additional pdf file supplementaryFile.pdf includes mathematical derivations and further descriptions of the models, explanation of methods used, additional results and extra figures.