Global attractors and extinction dynamics of cyclically competing species

Transitions to absorbing states are of fundamental importance in nonequilibrium physics as well as ecology. In ecology, absorbing states correspond to the extinction of species. We here study the spatial population dynamics of three cyclically interacting species. The interaction scheme comprises both direct competition between species as in the cyclic Lotka-Volterra model, and separated selection and reproduction processes as in the May-Leonard model. We show that the dynamic processes leading to the transient maintenance of biodiversity are closely linked to attractors of the nonlinear dynamics for the overall species’ concentrations. The characteristics of these global attractors change qualitatively at certain threshold values of the mobility and depend on the relative strength of the different types of competition between species. They give information about the scaling of extinction times with the system size and thereby the stability of biodiversity. We define an effective free energy as the negative logarithm of the probability to find the system in a specific global state before reaching one of the absorbing states. The global attractors then correspond to minima of this effective energy landscape and determine the most probable values for the species’ global concentrations. As in equilibrium thermodynamics, qualitative changes in the effective free energy landscape indicate and characterize the underlying nonequilibrium phase transitions. We provide the complete phase diagrams for the population dynamics and give a comprehensive analysis of the spatio-temporal dynamics and routes to extinction in the respective phases.

Figure 1

Typical spatial patterns in the May-Leonard model. Color (gray scale) denotes the concentration of species $A$, $B$, and $C$, with red (medium gray) signifying a site dominated by species $A$, green (light gray) a site dominated by species $B$, and blue (dark gray) being a site dominated by species $C$. (a) For large diffusion constants ($D=1.5\times {10}^{-3}$), the dynamics shows global oscillations with periodic switching between states with mainly one species present. (b) For intermediate values of the diffusion constant ($D=5\times {10}^{-4}$), we observe planar traveling waves. Here two of the domains take a characteristic domain size dictated by the diffusion constant. The third domain then occupies the rest of the system. (c), (d) For even smaller diffusion constants ($D=3\times {10}^{-4}$ and $D=6\times {10}^{-5}$), pairs of rotating spirals appear. The vertices of these spirals move very slowly on a time scale much larger than the time scale corresponding to the propagation speed of their arms. The wave length of the spirals decreases when the diffusion constant is reduced. The system size for all snapshots was $L=60$, and the carrying capacity for each site was chosen to be $M=8$.

Figure 2

Free energy landscapes of the May-Leonard model for various values of the diffusion constant $D$. The dynamics of the overall densities $\overline{a}$, $\overline{b}$, $\overline{c}$ is strongly confined to the invariant manifold of the well-mixed model, Eq. (17). To study the mechanisms underlying the distinct spatio-temporal patterns found in the spatial, stochastic May-Leonard model we projected the probability for the overall densities onto the invariant manifold of the rate equations (17) for different values of the diffusion constant $D$. Color (gray scale) denotes the logarithmic probability to find the system globally in a specific state before reaching one of the absorbing states, such that red (medium gray) denotes a high probability, yellow (light gray) a medium probability, and blue (dark gray) a low probability. The absorbing states themselves are not part of the statistics. For large $D$, no stable spatial structures can form and the dynamics corresponds to the well-mixed case, Eqs. (17). As $D$ becomes smaller than $D\approx 9.5\times {10}^{-4}$ an attractor of the global dynamics emerges, effectively stabilizing the system against extinction. This attractor corresponds to planar, traveling waves with oscillating overall densities, and grows in radius with decreasing $D$ due to a decreasing wavelength of the planar waves. For $D\le 3\times {10}^{-4}$, a second attractor emerges, corresponding to rotating spirals. As a result of a decreasing wavelength of the spiral patterns the second attractor's radius decreases with the diffusion constant, while the attractor corresponding to the traveling waves diminishes. Parameters were $L=60$ and $M=8$. Each plot was averaged over at least 100 realizations of the stochastic spatial dynamics.

Figure 3

Free energy landscape of the global dynamics of the two-dimensional May-Leonard system. The value of the Lyapunov function $\mathcal{L}$ is a measure of how close a specific state is to the boundary of the invariant manifold ($\mathcal{L}=0$). Color (gray scale) denotes the logarithmic probability to find the system at a specific value of $\mathcal{L}$, i.e., the free energy formally defined in Eq. (19). Red (medium gray) signifies small values of the free energy (minima of the potential) and thereby an attractor of the global dynamics. Yellow (light gray) denotes intermediate values, and blue large values of the free energy. The free energy landscape changes qualitatively at two threshold values for the diffusion constant $D$. For large values of $D$ the effective free energy has minima in the center ($\mathcal{L}=0.037$), where the dynamics starts, and at the boundaries of the invariant manifold ($\mathcal{L}=0$). At a first threshold $D_1\approx 9\times {10}^{-4}$ an attractor emerges, which moves away from the reactive fixed point ($\mathcal{L}=0.037$) with decreasing values of $D$. Below a second threshold, $D\approx 4.5\times {10}^{-4}$, a second attractor emerges near the reactive fixed point, coexisting with the first one. For even smaller mobilities the dynamics is solely determined by the attractor near the reactive fixed point. Comparing with our simulations we find that these attractors correspond to global oscillations (heteroclinic orbits), planar waves, and rotating spirals, respectively. The stochastic simulations were performed on a square lattice of linear size $L=60$ and with a carrying capacity $M=8$ for each site. For each values of $D$, the histogram was averaged over at least 100 realizations.

Figure 4

Mean lifetimes (dark dots) and coefficient of variation (gray squares) of the two-dimensional May-Leonard system as a function of the diffusion constant $D$. The mean lifetime increases abruptly at a first threshold value $D_1$ (indicated by a dashed line) where planar waves form. After passing through a maximum the lifetime decreases again. This is due to planar waves which become unstable with decreasing correlation length. At the second threshold $D_2$ (dashed line) rotating spirals become possible. With further decreasing diffusion constant the mean lifetime is asymptotically described by a power law $D^{-2}$. The coefficient of variation is a dimensionless measure for the dispersion of the probability distribution of $T$. The dispersion near the upper threshold $D_1$ becomes large, i.e., we observe dynamical regimes on a variety of different time scales.

Figure 5

Illustration of the instability of wave fronts in the cyclic Lotka-Volterra system. The initial condition at $t=0$ was chosen as three domains of equal size in cyclic order. The pictures show snapshots at times $t=135$, 180, and 300. Color (gray scale) denotes species concentrations as described in Fig. 1. Parameters were $D={10}^{-4}$, $M=8$, and $L=80$.

Figure 6

(a) Probability that the system with Lotka-Volterra reactions reaches an absorbing state before $t=N$. We observe a sharp transition from survival to extinction at $D\approx 3\times {10}^{-3}$. In the well-mixed limit the extinction probability converges to a finite value of 0.8 ($M=8$, $L=60$). (b) and (c) The scaling of the mean first passage time into any of the absorbing states with the system size $N$. Two phases can be identified. For $D<3\times {10}^{-3}$ the scaling becomes exponential, hinting at an escape process from a metastable state. In the well-mixed case ($D=100$) the scaling is linear, in agreement with the escape out of a neutrally stable state. In the intermediary regime our results are in agreement with both a logarithmic and a linear scaling.

Figure 7

Probability to find the system globally in a specific state before reaching one of the absorbing states. Red (medium gray) denotes a high probability, yellow (light gray) a medium probability, and blue (dark gray) a low probability. One observes the emergence of an attracting limit cycle of the global dynamics. The attractor grows in radius with increasing $D$ and eventually reaches the boundaries of the invariant manifold. For even larger values of the diffusion constant the attractor lies outside of the invariant manifold and the global dynamics is essentially neutral. The histograms were sampled over 100 trajectories until $T=N$. Parameters were $M=8$ and $L=60$.

Figure 8

The global phase portrait of the Lotka-Volterra system. For each diffusion constant $D$ we plot (in color code or gray scale) the probability to find the system before $t=N$ at a specific value of the Lyapunov function $\mathcal{L}=\overline{a}\overline{b}\overline{c}$. Red (medium gray) denotes a high probability, yellow (light gray) a medium probability, and blue (dark gray) a low probability. Three dynamical regimes corresponding to neutrally stable orbits, system-wide oscillations, and convectively unstable spirals can be identified and linked to their corresponding attractors of the global dynamics. Note that the attractor vanishes at $D\approx 3\times {10}^{-3}$ corresponding to the abrupt increase of extinction probabilities in Fig. 6. The simulations were done with $M$ $=$ 8 and $L=60$, and, for numerical reasons, stopped at $T=N$. For each of roughly 20 data points in $D$ we averaged over approximately 100 trajectories.

Figure 9

(a)–(c) Snapshots of the spatial distribution of species for different values of the fraction of Lotka-Volterra reactions $s$ indicated in the graph. Color denotes (gray scale) the concentration of the species $A$, $B$, and $C$, as described in Fig. 1. With increasing $s$ spirals become convectively unstable, i.e., they move, annihilate, and then appear again. (d)–(f) To illustrate the destabilization of spiral waves we computed for each lattice site the distance from the reactive fixed point $\left|\mathbf{y}\right(a,b,c\left)\right|$. Dark points show sites where $\left|\mathbf{y}\right(a,b,c\left)\right|$ is below a certain threshold, thereby indicating the position of spiral vertices. Parameters were $D={10}^{-4}$ (corresponding to the regime, where spirals and waves are possible in the May-Leonard model), $M=8$, and $L=80$.

Figure 10

Average first passage time into any of the absorbing states as a function of the diffusion constant $D$ and the fraction of Lotka-Volterra reactions $s$. Red (medium gray) denotes a large lifetime, yellow (light gray) a medium lifetime, and blue (dark gray) a small lifetime. For $s=0$ we obtain the mean lifetimes shown in Fig. 4. The dynamics is essentially governed by heteroclinic orbits, traveling waves, and rotating spirals. With increasing $s$ the planar waves become more and more stable and dominate the dynamics for a full order of magnitude in $D$. The prominence of traveling waves leads to a local maximum in the mean lifetimes. For even higher $s$ the system undergoes an Eckhaus instability (the analytical value is denoted by a dashed line), where planar waves become unstable. The dynamics is roughly comparable to the heteroclinic orbits in the May-Leonard model. Neutral orbits are driven to the boundary of the invariant manifold by a limited fraction of May-Leonard reactions. For $s=1$ we again recover the dynamics of the Lotka-Volterra model studied in Sec. 3. For each of the approximately 400 data points averages were taken over about 100 trajectories. Due to numerical constrains simulations were stopped at $T={10}^7$. Parameters were $M=8$ and $L=60$.

Figure 11

Probability to find the system in a specific state before reaching one of the absorbing fixed points for fixed $D=3\times {10}^{-4}$. The histogram is projected onto the invariant manifold of the rate equations (4). We varied the fraction of Lotka-Voltera reactions from $s=0$ (top left) to $s=1$ (bottom right). Starting at the classic May-Leonard model ($s=0$), where attractors for planar waves and rotating spirals can be identified, the attractor for the spirals disappears with growing $s$. The remaining attractor contracts to the reactive fixed point and also, when the system undergoes an Eckhaus instability, dissolves. For an even larger fraction of Lotka-Volterra reactions the global dynamics is driven outward by a limited fraction of May-Leonard reactions and is comparable to the heteroclinic orbits found in the May-Leonard model. When a majority of the reactions are of Lotka-Volterra type, i.e., for $s$ not much smaller than 1, we again observe the emergence of an attracting limit cycle corresponding to system-wide oscillations. Parameters where $M=8$ and $L=60$.

Figure 12

To study the most intriguing line of Fig. 10, $D=3\times {10}^{-4}$, in more detail we computed the effective free energy $\mathcal{F}$ at specific values of the Lyapunov function $\mathcal{L}$ for different values of $s$. Red (dark gray) denotes minima of the potential and thereby attractors of the global dynamics. Yellow (light gray) signifies intermediate values, and blue (dark gray) large values of the effective free energy. We identify several regimes depending on the relative strength $s$ of the different types of competition. For $s=0$ we recover the coexisting wave and spiral attractors of the May-Leonard model. With increasing values of $s$ only the wave attractor remains and approaches the reactive fixed point of the global dynamics ($\mathcal{L}=0.037$). At the Eckhaus instability (dashed line) the wave attractor dissolves. Instead, an attractor corresponding to global heteroclinic orbits emerges. Only when almost all reactions are of Lotka-Volterra type does an attractor close to the reactive fixed point emerge. The latter corresponds to the limit cycle found in the cyclic Lotka-Volterra model; see Sec. 3. Comparing with our simulations we find that these attractors are linked to rotating spirals, planar waves, global, heteroclinic orbits, and system-wide oscillations. Simulation parameters were $M=8$ and $L=60$.