Frequency-dependent fitness induces multistability in coevolutionary dynamics

Abstract

Evolution is simultaneously driven by a number of processes such as mutation, competition and random sampling. Understanding which of these processes is dominating the collective evolutionary dynamics in dependence on system properties is a fundamental aim of theoretical research. Recent works quantitatively studied coevolutionary dynamics of competing species with a focus on linearly frequency-dependent interactions, derived from a game-theoretic viewpoint. However, several aspects of evolutionary dynamics, e.g. limited resources, may induce effectively nonlinear frequency dependencies. Here we study the impact of nonlinear frequency dependence on evolutionary dynamics in a model class that covers linear frequency dependence as a special case. We focus on the simplest non-trivial setting of two genotypes and analyse the co-action of nonlinear frequency dependence with asymmetric mutation rates. We find that their co-action may induce novel metastable states as well as stochastic switching dynamics between them. Our results reveal how the different mechanisms of mutation, selection and genetic drift contribute to the dynamics and the emergence of metastable states, suggesting that multistability is a generic feature in systems with frequency-dependent fitness.

Figure 1.

The stationary distribution for the model system exhibits 3 maxima corresponding to metastable points. (a) Theoretical curve from equation (3.15) (solid line) in perfect agreement with data from simulations with N = 1000 (cross symbols). (b) Fitness functions of genotype A (solid line) and B (dashed line). Interaction parameters are , , and (see equation (3.1) and (4.1)) while the mutation rate is (). (Online version in colour.)

Figure 2.

Asymmetric mutation rates cause the emergence of a new maximum of the stationary distribution. (a) Stationary distribution of the Fokker–Planck equation for (solid line) and (dashed line) for a system with selective advantage for genotype B as shown in (b) (solid line: genotype A; dashed line: genotype B). Solid and dashed lines in (a) show the theoretical curves from equation (3.15), crosses the data from simulations with N = 1000 (cross symbols, ; plus symbols, ). Interaction parameters are , , and (see equation (4.1)) while mean mutation rate is . (Online version in colour.)

Figure 3.

The mutation rates strongly influence the system's dynamics. (a) Stationary solutions of an example system with interaction functions as in equation (4.1) for three different symmetric () mutation rates: (solid line), (long dashed line) and (short dashed line). This illustrates that for low mutation rates the dynamics stay close to the system's edges for long times. Here, differences in the mutation rates only slightly affect the shape of the stationary distribution. On the other hand, (b) shows that for high mutation rates () the dynamics are drawn stronger towards the middle. Here, differences in the mutation rates have a stronger impact, as the three curves with (solid line), (long dashed line) and (short dashed line) demonstrate. The increasing shifts the existing stable state towards the edge of the system until it vanishes. System parameters were N = 1000 and a symmetric interaction function according to equation (4.1) with one stable state at x = 0.5 with , . (Online version in colour.)

Figure 4.

The theoretically obtained stationary distribution well fits the data from simulations for large population sizes N. (a) Theoretical curve from equation (3.15) (solid line) for the same system as in figure 1. Data points show empirical distributions as defined in equation (4.2) obtained in simulations for N = 50 (cross symbols), N = 100 (asterisks), N = 500 (filled circles) and N = 1000 (plus symbols). (b) Distances and of simulation data and theoretical curve for different N for the mean distance measure defined in equation (4.3) (cross symbols) and the maximum distance measure defined in equation (4.4) (asterisks), demonstrating that the distance decays for increasing N. The solid line and the dashed line are added as a guide to the eye for the relation between measured distances and system size N. The measured empirical distributions for both (a,b) were obtained by simulating the system dynamics from an initial state drawn from for a mixing time and then recording the density for a time . (Online version in colour.)

Figure 5.

The stationary distribution (3.15) obtained from the Fokker–Planck equation well fits the exact solution (4.6) from the Master equation for weak selection, but not for strong selection. (a) Stationary solution from the Fokker–Planck equation (solid line) for the system with interaction functions defined by equation (4.5) for very weak selection together with the solution from the Master equation (cross symbols). The distributions of the same system for weak selection and strong selection are shown in (b,c), respectively. All curves are plotted logarithmically to better allow a comparison of the deviations over all scales. (c) The curve from the Fokker–Planck equation does not fit the exact Master solution very well. This is quantified in (d) showing the mean distance between Fokker–Planck and Master solution (cross symbols) and between Fokker–Planck solution and an empirical density obtained from simulations (asterisks). The system parameters were N = 1000, , and . The measured empirical distributions were obtained by simulating the system dynamics from an initial state drawn from for a mixing time and then recording the density for a time . (Online version in colour.)

Figure 6.

Multiple stability dynamics for periodic interaction functions. (a) Theoretically computed stationary distribution from equation (3.15) (solid line) together with simulation data with N = 1000 (cross symbols) for the interaction functions given in (4.7). The computed stationary distribution diverges for and while the simulation data remain finite owing to the finite number of individuals. (b) Sample path that exhibits switching between the different maxima of the distribution. Parameters are , , and . (Online version in colour.)