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.


S1 Additional results
Selection intensity: Here, we analyse the impact of the selection intensity w which modulates the pulling force in models with an attractive fixed point. For a more robust result we compare the simulations (lines in Figure S1) with sojourn times (section S6) calculated for the discrete time dtEvo + and dtEvo processes. Since only discrete time and discrete state processes can be represented by a transition matrix, we can only apply this approach to the discrete time constant population size processes. Although the Evo + process has an attracting fixed point, intuitively making extinction times longer than in the Evo process, for low population sizes and weak selection we see the opposite -fast extinction. Furthermore, while the extinction time increases exponentially with increasing intensity of selection in the Evo + process, the Evo extinction times stay comparably constant, which is not surprising considering the neutral stability. One interesting and unexpected result is that extinction time is lower for increasing selection intensity of one species while keeping the other constant. This occurs when one of the selection intensities is very low: some dark lines are decreasing, especially for example w H = 0.2, and colours are reversed for low values of w P .
Dimensionality: We next include more types to observe the exponential decline of diversity, here simply the number of genotypes n H and n P alive ( Figure S2). The Evo + process simulated with a Gillespie algorithm is updated such that the interaction matrix is normalised depending on the number of strains (otherwise there is an imbalance between matching and non-matching pairs which results in change of selection strength). Oscillating selection can give an advantage to any type, but at different time points. A previously extinct parasite type (P 2 ) is reintroduced manually at time point 16, where the abundance of the corresponding host is especially high.   Figure S2: Diversity decline of subtypes of hosts and parasites. Example of an Evo + process implemented with a Gillespie algorithm. The simulations start with equal abundance of all 20 types H i (0) = N H /20 and P i (0) = N P /20. At time point t = 16 a well adapted but extinct P 2 = 1 is reintroduced manually, while P 5 is reduced by one individual to keep N P constant. Parameters: N H = 200, N P = 200, w H = w P = 1, α = 1, β = 0. The legend shows host and parasite subtypes present at t = 16.

S2 Model descriptions
Verbal model descriptions and some parameter explanations are available in the main text, here we show the mathematical details.

S2.1 The Evo + and Evo process
In the Evo + process (often Moran process in evolutionary game theory) and Evo process (often pairwise comparison or local update process) population size N H and N P is constant, so it is enough to focus on one type H = H 1 = N H − H 2 and P = P 1 = P N − P 2 . Then the fitness of the second type has a negative influence on the first type and we can write fitness as f s r , with s ∈ {+, −} indicating whether the fitness has a positive or negative influence on the focal type and r ∈ {h, p} denoting the population. Fitness depends on the payoff taken from the game played between hosts and parasites with relative abundances h := H/N H and p := P/N P . The infection matrices M p = α β β α The strength of the influence of this payoff on fitness is modulated by the selection intensity w r ∈ [0, 1]. Here we use a linear dependence In an Evo + process the local probabilities depend on the fitness of a type divided by an average fitness (fitness of subtypes weighed with the relative abundance) of the species. (2) For the Evo process, the local transition probabilities are influenced by the fitness difference between the two types scaled with the maximal possible payoff difference where ¬s is the other type (chosen for death). We often use α = 1 and β = 0, then ∆π rmax = 1, else ∆π rmax = α − β. The local rates Φ s r are the specific rates for the birth-death reactions. Multiplied with the relative abundance of the reactants, we get the global transition probabilities T s r = r (1 − r) Φ s r . In the discrete time processes host and parasite populations are updated simultaneously and all transition probabilities sum up to one. For each population, there must therefore be a transition probability for increasing, decreasing or not changing the number of the focal type. The respective probabilities are T + r , T − r and T 0 r = 1−T + r −T − r . For the two interacting populations this yields nine probabilities, Tab. S1. In the computer simulation, in each time step ∆t = 1, a random number is compared with the cumulative probabilities and the resulting transition is carried out by adjusting the population accordingly and updating time to t + 1.

S2.2 Gillespie algorithm
In a Gillespie algorithm we assume independent interactions with one or two reactants that react with a certain rate, which is the product of the reactant concentrations multiplied with a reaction constant (not Table S1: Transition probabilities T l (h, p) necessarily a constant here). If reactions happen with waiting times that follow an exponential distribution we can calculate the time of the next reaction, for each reaction with mass-action rate φ and the help of a random number R ∈ [0, 1], We then choose the reaction that happens first and update population and time.
For reactions with one reactant, the reaction rates are multiplied with the absolute abundance (e.g. φ = H 1 b for the first reaction in the EcoEvo + process), but when two reactants are necessary, the probability of that reaction would be amplified. The reaction is therefore scaled. In the Evo + /Evo process pair reactions are scaled with N H (N P ) for the first two (last two) reactions (see Table S2), so that φ = H 1 H2 N H Φ + H , etc. In the EcoEvo + and EcoEvo process reactions with two reactants are scaled with the carrying capacity K (φ = H 1 H2 K µ for the third reaction). The last case is very close to the usual approach in chemistry, where reactions rates with more than one reactant are scaled with the system volume.

S2.3 EcoEvo + and EcoEvo process
For the completely independent EcoEvo models reactions with two reactants are scaled with an extrinsic carrying capacity K. These are reactions where host and parasite interact with rate λ = λ0 K or where hosts are in competition with each other µ = b h K (or µ = 0 in the EcoEvo model). The choice of µ becomes apparent in the deterministic limit, where the carrying capacity in the logistic term is K = b h µ (see Tab. S3). Note that the carrying capacity does not equal the population size of the hosts when parasites are present. The initial conditions are chosen to be similar with the game theory processes by varying the value of the fixed point through b h .

S2.4 Hybrid process
Like the Evo + an Evo process the Hybrid model takes its interaction rates via a game theory approach from the payoff matrix. In the Evo + and Evo process the payoff matrix for the host is adjusted to realise positive fitness. Here, the set-up is closer to typical host-parasite modelling by using only one payoff matrix for positive rates b P1 , b P2 for the parasite and negative rates d H1 , d H2 for the host. With M = α β β α . When w H = w P = 1 the model is equivalent to the independent reactions with two reactants, however, scaled with different quantities H 1 + H 2 or P 1 + P 2 , which are dynamic. This model is in between the EcoEvo and the Evo + /Evo processes.
death host and birth parasite when w H = w P = 1 : S3 Stochastic differential equation and deterministic limit S3.1 Master equation In the limit of N H , N P → ∞ and ∆t → 0 the stochastic process is described by a Fokker Planck equation with a selection term (the deterministic ordinary differential equation, often called drift in physics) and a diffusion term (stochastic term often called genetic drift in biology). The Fokker Planck equation can also be converted into a stochastic differential equation. As an example we explain this procedure for the discrete time Evo + or Evo process.
From the transition probabilities above, one can write down the master equation The probability P (h, p, t + ∆t) of being at the specific state (h, p) at time t + ∆t is the probability of being in that state plus the probability of being in a state close by at time t multiplied with the probability of transitioning to the new state (h, p) minus the probability of leaving that state. See Tab. S1.

S3.2 Fokker-Planck equation
We can now use the shift operator E ∆x on the left hand side and on the right hand side of the master equation. Then approximating to first order and to second order ( 1 N 2 , see (van Kampen, 1997;Risken, 1996;Gardiner, 1985)) we get The resulting Fokker-Planck equation emerges when we collect similar terms The drift vector a = (a 1 , a 2 ) and the diffusion matrix D = d1 d1,2 d2,1 d2 are calculated as follows, where we have dropped the arguments (h, p) for readability.
We have summarised the results for the infection matrices M p = α β β α and M h = β α α β in Tab. S3.

For the Evo
1−w P +w P π P and Φ + P + Φ − P = 2−2w P +w P (α+β) 1−w P +w P π P . And for the Evo process Φ + We have normalised the reaction rates in the discrete time processes so that all probabilities sum up to one. The other processes' reaction rates are on the order of N H or N P (also see appendix section S2.2), for the continuous time processes the Fokker Planck equation is therefore faster. Furthermore, some of the processes are four dimensional with drift vector a = (a 1 , a 2 , a 3 , a 4 ) for H 1 , H 2 , P 1 , P 2 and diffusion matrix D = d1 0 0 0 0 d2 0 0 0 0 d3 0 0 0 0 d4 . Here, all correlative noise (non-diagonal entries) is zero because all the reactions only change one reactant at a time. Note that without diffusion, the Hybrid process is equivalent to the scaled two-dimensional replicator dynamics because a 1 + a 2 =Ḣ 1 +Ḣ 2 = 0 which implies H 1 + H 2 = N H .

S3.3 Stochastic differential equation
From the Fokker Planck equation a stochastic differential equation (SDE) can be constructed. The SDE is a differential equation with a deterministic part, which has to be integrated with respect to time and a second term which is integrated with the Wiener process (continuous time random walk/Brownian motion). While the Fokker Planck equation describes the probabilities for certain states in time, and thus a probability distribution, the SDE accounts for the change of a random variable. Integrated with respect to one particular realisation of a random walk (Wiener process), it will give one particular realisation of the stochastic process.
To derive a non-unique SDE from the Fokker Planck equation we take the "square root" of the matrix D = d1 d1,2 d2,1 d2 = B T B with d 1,2 = d 2,1 . A Cholesky decomposition can achieve this, where c(x) is the complex conjugate of x. This is only important if d 1 can become negative which is not the case for our purposes. With the drift term a(X t ) diffusion matrix B = B(X t ) we set up the stochastic differential equation in the Itô sense where the stochastic variable X t = h p describes the relative abundance of the first host and first parasite, and W is the Wiener process. The terms a 1 and a 2 are equivalent to the differential equation (deterministic terms in Table S3) and describe the N → ∞ limit of the stochastic process, which is e.g. the replicator equation for the Evo process and the adjusted replicator dynamics for the Evo + process.
The Fokker Planck equation is a deterministic partial differential equation, which can be solved in simple cases. Although this is not the case here, we can learn from the matrix D how the variance in the distribution is affected by the parameters. Again, although all seven models stem from the same biological mechanisms, the results are quite different. The Evo + , Evo, EcoEvo + and EcoEvo processes have only diagonal elements, therefore the random term is not correlated across species (2D models) or even types (4D models).
Stochastic differential equations can be numerically integrated in a robust way (Rößler, 2010;Kloeden and Platen, 1992) using random numbers for the Wiener process. The advantage of an SDE with respect to a stochastic simulation (via Gillespie's algorithm) is the independence of run time on the population size N . In a stochastic simulation each individual in the population is simulated. An SDE, like an ODE is integrated stepwise, thus the limiting factor here is the time step used. For small population numbers as displayed in the figures of the main text both algorithms are comparable in speed. As population size increases stochastic simulations become infeasible, yet SDE integration does not increase the runtime. Model deterministic terms stochastic terms

S4 Deterministic properties
Fixed points occur when the dynamics have come to a halt and there is no more change in (relative) abundances. In the two-dimensional deterministic models (dtEvo + , Evo + , dtEvo, Evo and Hybrid -if normalised), which can be described by the replicator dynamics or the adjusted replicator dynamics, there are trivial points, where one subtype is extinct: (h = 0 or h = 1) and (p = 0 or p = 1). The less trivial fixed point is the inner coexistence fixed point with h * = p * = 0.5 in the symmetric case where there is no distinction between subtypes except the specificity to a matching host/parasite. The fourdimensional deterministic models (EcoEvo + and EcoEvo) have coexistence fixed points H * i = dp λ = dp K λ0 there is no logistic growth. A linear stability analysis is carried out to learn more about the properties of the fixed point. The Jacobian at this inner fixed point J(h * , p * ) = dḣ dh dḣ dp dṗ dh dṗ dp h=h * ,p=p * is calculated to reveal the characteristic function det (J(h * , p * ) − λI) = 0. If the real parts of all Eigenvalues are negative, the fixed point is attractive and thus (in the terms of evolutionary game theory) an evolutionary stable strategy ESS. If at least one Eigenvalue has a positive real part, the state is a repellor (unstable). And if all real parts are zero, it is not enough to conduct a linear stability analysis. For the replicator dynamics and the adjusted replicator dynamics we find only imaginary Eigenvalues where the real part is zero, indicating neutral stability. This can be addressed by more sophisticated methods.
A constant of motion (Hofbauer, 1996;Koth and Sigmund, 1987) for the replicator dynamics reveals that the dynamics are volume preserving. The fixed point is neutrally stable. All trajectories follow orbits as shown in Fig. S3. Also, from Hofbauer and Sigmund (1998, chapter 11) we know that in our example with M p = α β β α and M h = β α α β we have a c-zero-sum game. Then the fixed point is globally asymptotically stable (attractive) in the adjusted replicator dynamics. A simple stability analysis shows that the EcoEvo + is attractive for many parameters but neutral when µ = 0 in the EcoEvo interactions with no constraint on population size.

S5 Approximate diffusion for a constant of motion
For the neutrally stable models that lead to replicator dynamics (dtEvo, Evo, Hybrid) and to the original independent reactions dynamics (EcoEvo) one can set up a constant of motion H. Here we use the replicator dynamics as an example.
A constant of motion is a function of the state variables in the system that does not change in time dH dt = 0 when a trajectory is not perturbed by chance. Intuitively one can think of a system with constant energy where the dynamics are confined to the constant energy contour lines. For the two-dimensional replicator dynamics and a matching allele infection matrix M p = ( 1 0 0 1 ) and M h = ( 0 1 1 0 ) the constant of motion is h, p ∈ [0, 1] are the relative abundances of host type one and parasite type one (1 − h and 1 − p are the relative abundances of type two) and the constant of motion is scaled so that it has the value one when both types are equally abundant (h = p = 1/2) and the value zero when any type is extinct (Claussen, 2016). The deterministic trajectory would start with a certain value of H, depending on the initial condition, and would follow that line indefinitely (see Fig. S3). In a stochastic model the constant of motion loses this meaning, but it can still be used as an observable to measure the distance of the state to the inner fixed point or the outer border where a type goes extinct. Claussen (2007) uses the average change of this observable to analyse under what conditions the system can be forced to return to the inner fixed point (drift reversal) or whether extinction is inevitable owing to the outward diffusion in a "battle of the sexes" game. In fact, the quantity can be used as an observable also for models that are not neutrally stable, but where H = 0 still means extinction and h, p are the two dimensions. Hence, we can use the following approach not only for the Evo process (which is equivalent to the replicator dynamics in a limit) but also the Evo + process (which is equivalent to the attractive adjusted replicator dynamics). The discrete time versions of these processes allow us to focus on the probabilities only, without having to keep track of time, since ∆t = 1. We now use the average change of this observable to analyse whether the system can be forced to return to the inner fixed point or whether extinction is inevitable due to the outward diffusion.
When a transition occurs, the value of the observable H also changes. For each transition T l (h, p) with l = 1, 2, ..., 9 and corresponding ∆h and ∆p there is a ∆H l (h, p) := H(h + ∆h, p + ∆p) − H(h, p). By weighing all possible changes in H with the corresponding transition probabilities one receives the average change of the observable for each configuration of host and parasite populations.
This is the average change in H for each coordinate (h, p) in each time step, but the actual extinction time remains to be calculated. The average diffusion over the whole space would be the integration over h, p ∈ [0, 1], which is not possible in all generality. It becomes possible when taking the weak selection approximation via a Taylor approximation around w H = w P = 0 is up to second order in w i , i ∈ {H, P }. After this approximation it is possible to integrate the formula to get a (phase space) average of the average diffusion.
For the matching allele infection matrices M H = ( 0 1 1 0 ) and M P = ( 1 0 0 1 ) in a discrete time Evo or Evo + process the average diffusion is with per time step. Now we can find out the time steps it takes to reach a certain value of H f inish starting in a certain value of H start .
We start in H start = 1, which is the coexistence fixed point with a maximum variability of host and parasite subtypes. The value of H when one subtype dies out is H f inish = 0.
This method only works if the average diffusion is negative. For example, with w H = w P = 1 and N H = 250, ∆H dtEvo is always negative but ∆H dtEvo + becomes positive for N P ≥ 85. The stronger the selection intensity in the dtEvo + process, the more the pull towards the coexistence (attractive) fixed point. For the dtEvo, the selection intensity has hardly any influence (there are only terms in 1 N 2 i , whereas in the dtEvo + there are terms in 1 Ni ). One has to be careful to use this formula since it was derived in a weak selection limit and large N . Yet the simulations present results for small N and relatively large w. However, we see that the order of magnitude is the same as in the simulated extinction times above. Yet the logarithmic scale, although convenient to represent the different orders of magnitude, may be misleading.

S6 Exact sojourn times
The constant population size processes with discrete time (dtEvo + and dtEvo) allow for a description of transitions as probabilities in each discrete time step ∆t = 1. For each state with index k ∈ {1, 2, ...(N H + 1)(N P + 1)} in the two-dimensional state space H 1 ∈ {0, 1, ..., N H } × P 1 ∈ {0, 1, ..., N P } there is a probability of transitioning to every other state (most of these probabilities are zero, because only one step is taken in state space for each time step, see also Table S1). This transition matrix Π ∈ [0, 1] (N H +1)(N P +1) can be ordered into a block matrix containing one block for transient states Q, one block for transitions from transient states to absorbing states R, and one identity matrix block for absorbing states I. Then the transition matrix is and the recursion equation for the probability density ρ(t) = (ρ k ) k over all states k ∈ {1, 2, .., (N H + 1)(N P + 1)}. Grinstead and Snell (2012) proved that the time spent in the inner states, which are part of the matrix Q, depending on the initial state is an entry in (Q − I) −1 . Hindersin et al. (2016) apply this by solving the equivalent system of equations (Q − I)x = −1 for x. We have applied this numerical approach to the inner states H 1 ∈ {1, 2, ..., N H − 1} × P 1 ∈ {1, 2, ..., N P − 1} which then gives the exact average time to extinction of one of the subtypes in one of the populations.

S7 Literature overview
In the main text, Table 1 summarises some mathematical models on host-parasite dynamics found in the literature. We have based the choice of our models on the visibility of the studies, the relatedness and comparability to our work and the availability of the information found in the publications, that is necessary to complete the table. We have aimed at focussing on explicit host-parasite models (HP), that build on the classic Lotka-Volterra or Rosenzweig-MacArthur models or the evolutionary game theory literature. There are many such publications from before the year 2000 (Schaffer and Rosenzweig, 1978;Seger, 1988;Nee, 1989, and more), but we aimed at comparing predominantly the more recent literature. However, many models are host-focussed, such as susceptible-infected (SI) epidemiological models. Many host centred models include parasite genotype frequencies and infection probabilities, but the host dynamics are more complex than the parasite dynamics, and the focus is often on host evolution or the maintenance of sexual reproduction in the host. We have briefly discussed these studies, yet we are predominantly interested in co-evolution models, where parasites play a similarly important role (HP models). Furthermore, we have also included some studies that are concerned with slightly less related topics to show the diversity in assumptions and authors that this field provides. Some of these models include stochasticity and a changing population size (spatial structure model of Boots and Sasaki (1999)) while others include stochasticity but not a variable population size (model on special host or parasite behaviour Abou Chakra et al. (2014)) and others include a variable population size but are deterministic (non-host-parasite Red Queen dynamics Bonachela et al. (2017)). This is, naturally, a subjective approach, and we have not excluded any publications for belligerent reasons. We also wish to caution, that this study is not a systematic literature review and we do not make the claim to have collected all relevant publication in this extensive field.