Death by starvation in May-Leonard models

TL;DR

This paper investigates how individual death due to starvation affects the dynamics and coexistence in spatial stochastic May-Leonard models with mutual predation, using 2D simulations with two and three species.

Contribution

It introduces the impact of starvation-induced death on population dynamics and coexistence in May-Leonard models, a novel aspect not previously explored.

Findings

Starvation death significantly influences population stability.

Coexistence is preserved through starvation mechanisms.

Simulation results highlight the importance of individual death thresholds.

Abstract

We consider the dynamics of spatial stochastic May-Leonard models with mutual predation interactions of equal strength between any two individuals of different species. Using two-dimensional simulations, with two and three pecies, we investigate the dynamical impact of the death of individuals after a given threshold number of successive unsuccessful predation attempts. We find that the death of these individuals can have a strong impact on the dynamics of population networks and provide a crucial contribution to the preservation of coexistence.

Figures

Equations

$$ij\to i\otimes ,$$

$$i\odot \to \odot i,$$

$$i\otimes \to ii,$$

$$F(k_W[l])=\sum _{k=1}^{N_x}{\left(S_{kl}-ST[k-k_W]\right)}^2$$

$$W(t)=\sqrt{\frac{1}{N_y}\sum _{l=1}^{N_y}{\left(k_W[l]-⟨k_W⟩\right)}^2},$$

$$⟨k_W⟩=\frac{1}{N_y}\sum _{l=1}^{N_y}{k}_W[l].$$

$$p=r(1-{\rho }_D),$$

$${\rho }_D=1-\frac{p}{r}.$$

Full text

1 Introduction

Competition is ubiquitous in nature, playing a fundamental role on the regulation of biodiversity. It is also a major driving force behind evolutionary change through natural selection. The simplest competition models, inspired in the pioneering work by Lotka and Volterra, and May and Leonard [1, 2, 3], consider the dynamics of two or three species subject to interspecific predation (or selection), mobility and reproduction interactions. Despite their simplicity, these models (see [4, 5] for recent reviews) incorporate some of the main ingredients responsible for the observed dynamics of many biological systems, and are able to reproduce some the dynamical features of specific biological systems with a limited number of species [6, 7, 8, 9].

More complex competition models, involving more species [10, 11, 12, 13, 14] and/or additional interactions [15, 16, 17, 18, 19], have also been investigated in recent years, revealing a plethora of complex dynamical spatial structures [20, 21, 22, 23, 24], diverse scaling regimes [14, 25] and phase transitions [26, 27, 28, 29, 30, 31, 32, 33, 34]. In some of these competition models species coexistence may last for an arbitrary amount of time, while in others it is transient. Coexistence-promoting mechanisms, responsible for maintaining coexistence over long periods of time, are usually associated to the ability of the species to increase (decrease) their population in response to negative (positive) perturbations to their typical abundances [35]. Among these, density-dependent mortality [36, 37] has been claimed to have a positive impact in promoting species coexistence (see also [38, 39, 40] for a discussion of the impact of density-independent mortality).

Here, we investigate a sub-class of spatial stochastic May-Leonard models characterized by mutual predation interactions of equal strength between any two individuals of different species. In their standard version, the dynamics of these models results in a network of one-species domains whose dynamics is curvature driven, with the characteristic size of the network of one-species domains growing proportionally to $t^{1/2}$ [13] — $t$ being the physical time. However, in practice, this growth is limited by the size of the simulation boxes, thus resulting in a limited period of coexistence.

The main aim of this letter is the determination of the impact on population dynamics of the death by starvation of individuals after a given number of successive unsuccessful predation attempts. We shall start by introducing our set of models in the following section. We then study the impact of death by starvation on the dynamics of initially flat and circular interfaces between spatial domains occupied by individuals of competing species, as well as the two-dimensional dynamics of population networks starting from random initial conditions. We will show that, under certain conditions, death by starvation prevents the endless growth of the characteristic length scale of the network of one-species spatial domains, acting as a coexistence-promoting mechanism.

2 Models

In this letter we shall investigate the dynamics of May-Leonard models with mutual predation interactions of equal strength between any two individuals of different species. To this end, we shall perform square lattice simulations with periodic boundary conditions in which each one of its ${\mathcal{N}}$ sites may be either empty or occupied by a single individual. The species are labelled by the number $i$ (or $j$), with $i,j=1,...,N$ — in this letter we shall only consider models with two or three species ($N=2$ or $N=3$). Empty sites shall be denoted by $\otimes$. The number of individuals of the species $i$ and the number of empty sites will be denoted by $I_{i}$ and $I_{\otimes}$, respectively — the density of individuals of the species $i$ and the density of empty sites shall be defined by $\rho_{i}=I_{i}/{\mathcal{N}}$ and $\rho_{\otimes}=I_{\otimes}/{\mathcal{N}}$, respectively. The possible interactions are predation

[TABLE]

mobility

[TABLE]

and reproduction

[TABLE]

where $\odot$ represents either an individual of any species or an empty space.

Mobility, reproduction and predation interactions occur with probabilities $m$, $r$ and $p$, respectively (the same for all species). For the sake of definiteness, we shall take $m=0.5$ and $m+p+r=1$ in all the simulations. At each simulation step, the algorithm randomly picks an occupied site to be the active one, randomly selects one of its four adjacent neighbour sites to be the passive one, and randomly chooses an interaction to be executed by the individual at the active position. These three steps are repeated until a possible interaction is selected. If predation is selected, the impossibility of executing an interaction happens when the passive is an empty site or the passive and active individuals are of the same species. On the other hand, if reproduction is selected, the interaction only takes place if the passive is an empty site. Whenever the selected interaction is not executed, the active individual is said to have carried out an unsuccessful interaction attempt. A generation time (our time unit) is defined as the time necessary for $\mathcal{N}$ successive (and successful) interactions to be completed.

The non-standard ingredient in our simulations is the death of an individual after a given number ${\mathcal{N}}_{u}$ of successive unsuccessful predation attempts — in our model the most recent number of successive unsuccessful predation attempts of its progenitor is passed on to every newborn individual at the time of birth. This means that the ability of a newborn to survive unsuccessful predation attempts is strongly dependent on the strength of its progenitor at the time of birth (the strength being defined as ${\mathcal{N}}_{u}$ minus the latest number of unsuccessful predation attempts).

3 Dynamics of initially flat/circular interfaces

In this section we study the impact that a finite value of ${\mathcal{N}}_{u}$ has on the dynamics of spatial stochastic two-species May-Leonard models with mutual interspecific predation interactions of equal strength. To this end, a large number of spatial stochastic numerical simulations has been performed with the following parameters: $m=0.5$, $r=0.3$, $p=0.2$, and ${\mathcal{N}}_{u}=25$ (any individual dies after 25 successive unsuccessful predation attempts) or ${\mathcal{N}}_{u}=\infty$ (standard case, no deaths by starvation). In this section we shall consider the dynamics of initially flat and circular interfaces, before studying the dynamics of population networks with random initial conditions in the following section.

3.1 Initially flat interface

Here we consider the dynamics of an initially flat interface separating the left and right halves of the lattice which are initially fully occupied by individuals of the red (1) and blue (2) species, respectively. Our simulations are performed on a $N_{x}\times N_{y}$ lattice. The position ($k_{W}(l)$) of the interface for each row $l$ may be found as the value of $k_{W}(l)$ that minimizes the sum

[TABLE]

for each value of $l$. Here, $ST[x]$ is the step function, defined as $ST=1$ for $x<0$, $ST=0$ for $x=0$, and $ST=-1$ for $x>0$, and $S_{kl}=0$, $S_{kl}=+1$, or $S_{kl}=-1$ depending on whether the site of coordinates $(k,l)$ is empty, occupied by an individual of the red (1) species, or occupied by an individual of the blue (2) species. We shall follow refs. [41, 42, 43], and define the interface width as

[TABLE]

where

[TABLE]

Figure 1 illustrates the evolution of the width of an initially flat interface obtained from $1000$ realizations considering ${\mathcal{N}}_{u}=\infty$ (upper panel) and ${\mathcal{N}}_{u}=25$ (bottom panel). Notice that, after an initial transient stage, the growth of the width of the interfaces for ${\mathcal{N}}_{u}=25$ is much faster compared to that of the standard ${\mathcal{N}}_{u}=\infty$ case — the evolution approaching a scaling regime with $W\propto t^{0.18}$ at large $t$ for ${\mathcal{N}}_{u}=\infty$. This much faster growth may also be confirmed in the snapshots, obtained for single realizations of the ${\mathcal{N}}_{u}=\infty$ and ${\mathcal{N}}_{u}=25$ models, shown in the inset panels. The inset panels show the development of complex dynamical structures along the interfaces for ${\mathcal{N}}_{u}=25$, leading to rougher interfaces compared to the standard ${\mathcal{N}}_{u}=\infty$ case. It is also clear from the inset panels that for ${\mathcal{N}}_{u}=25$ death by starvation results in a relatively low constant average density of individuals away from the interfaces — the average density being reached when the equilibrium between the mortality and reproduction rates is attained. This is the key property of models with finite ${\mathcal{N}}_{u}$, which is responsible for the growth of the interface thickness and roughness with time observed in Fig. 1.

Figure 2 shows the average evolution of the width of an initially flat interface with time obtained from $1000$ realizations with ${\mathcal{N}}_{u}=\infty$ (upper panel) and ${\mathcal{N}}_{u}=25$ (bottom panel), considering an alternative definition of interface width. In this case, the interface width $W_{*}$ defines the interval of $k$ for which the abundance of the two species is bellow $68\%$. Fig. 2 shows that, not only $W_{*}$ is much less sensitive to the box size than $W$, but also that the growth of $W_{*}$ with time observed in the ${\mathcal{N}}_{u}=25$ case is absent for ${\mathcal{N}}_{u}=\infty$. For ${\mathcal{N}}_{u}=25$ the evolution approaches a scaling regime with $W_{*}\propto t^{0.3}$ at large $t$.

Although the two definitions of interface width are physically distinct — $W_{*}$ being much less sensitive than $W$ to small wavelength fluctuations which do not introduce large modifications to the interface profile at each row — both evolve very differently in the ${\mathcal{N}}_{u}=25$ and ${\mathcal{N}}_{u}=\infty$ cases, thus capturing the impact of death by starvation on interface dynamics.

3.2 Initially circular interface

Let us now consider the dynamics of an initially circular interface. Figure 3 illustrates the evolution of the average density $\rho_{1}$ of individuals of the the red species ($1$) obtained using $1000$ realizations of the evolution of an initially circular domain containing individuals of that species surrounded by an outer domain containing individuals of the blue species ($2$), considering ${\mathcal{N}}_{u}=\infty$ (upper panel) and ${\mathcal{N}}_{u}=25$ (bottom panel). The solid red line shows the average value of $\rho_{1}$ as a function of the number of generations $t$ while the shaded region represents the sample standard deviation.

The upper panel of Fig. 3 shows that if ${\mathcal{N}}_{u}=\infty$ the initially circular domain always collapses. In this case, a standard curvature dominated regime is recovered, with the domain wall area being roughly proportional to $t-t_{c}$, for $t\leq t_{c}$ ($t_{c}$ being the time of collapse). This explains the approximately linear dependence of the average density of the inner species (1) on time.

The bottom panel of Fig. 3 shows that, contrary to what happens for ${\mathcal{N}}_{u}=\infty$, in the ${\mathcal{N}}_{u}=25$ case initially circular domains do not collapse despite the existence of an initial shrinking stage. Instead, the inset panels, displaying snapshots of a single realization taken at three different times, show that for ${\mathcal{N}}_{u}=25$ there is a significant departure of circular symmetry and that, after an initial shrinking stage, the initially circular domain may grow and split into separate subdomains.

4 Network Simulations

In this section we shall consider the results of spatial stochastic numerical simulations with random initial conditions in two spatial dimensions. At the beginning of the simulation each site is either occupied by an individual of any of the $N=2,3$ species or left empty with a uniform discrete probability of $1/(N+1)$ (except if stated otherwise, the following parameter values are assumed: $m=0.5$, $r=0.3$, $p=0.2$). The results of $1000^{2}$ simulations of the two and three species models for ${\mathcal{N}}_{u}=\infty$ (upper panels) and ${\mathcal{N}}_{u}=25$ (bottom panels) are illustrated in Figs. 4 and 5, respectively. The snapshots were taken after $1000$ (left panels), $4000$ (middle panels), and $16000$ generations (right panels).

In the absence of death by starvation (${\mathcal{N}}_{u}=\infty$) there are almost no empty sites deep inside the domains. Empty sites are created only when predation occurs on the interface between competing domains. Empty sites can move around as a result of mobility interactions, but are eventually filled as a result of reproduction. In the ${\mathcal{N}}_{u}=\infty$ case, the dynamics is curvature dominated, with the velocity of the interfaces being roughly proportional to their curvature. This has been shown to lead to a population network evolution whose characteristic lengthscale $L$ is roughly proportional to $t^{1/2}$ [13]. This growth of the characteristic scale of the network can be observed in the upper panels of Figs. 4 and 5, for simulations with two and three species respectively. Eventually, for ${\mathcal{N}}_{u}=\infty$ the size of the domains becomes of the order of the box size and the coexistence is lost (as illustrated in Fig. 7 considering a simulation on a smaller grid).

On the other hand, the bottom panels of Figs. 4 and 5 show that, in the presence of a mortality rate associated to insufficient predation (${\mathcal{N}}_{u}=25$), the average density of individuals has a peak in the interface regions, decreasing towards the interior of the domains, reaching an asymptotic value determined by the equilibrium between death and reproduction. This is responsible for the dynamical behaviour of initially flat and circular interfaces found in the previous sections and is the reason why the characteristic scale of the network does not change significantly from $t=1000$ to $t=16000$, as observed in the bottom panels of Figs. 4 and 5.

As discussed before, for ${\mathcal{N}}_{u}=\infty$ the empty sites are mainly concentrated on the borders of competing domains. The thickness of the interfaces of empty sites is roughly constant, which implies that the total interface length $L_{T}$ is roughly proportional to the number of empty spaces $I_{\otimes}$. On the other hand, the number of domains is roughly proportional the ratio between the total area $A$ and the average domain area $L^{2}$, where $L$ is the characteristic length scale of the network. It is also proportional to the ratio between the total interface length $L_{T}$ inside the simulation box and the average domain perimeter (which is proportional to $L$). Therefore, $A/L^{2}\propto L_{T}/L$, thus implying that $L\propto A/L_{T}\propto{\mathcal{N}}/I_{\otimes}=\rho_{\otimes}^{-1}$ [13]. Hence, for ${\mathcal{N}}_{u}=\infty$ the characteristic length scale of the network is inversely proportional to the number of empty spaces.

Figure 6 depicts the evolution of $\langle\rho_{\otimes}^{-1}\rangle$ (the angle brackets represent an ensemble average) with time for ${\mathcal{N}}_{u}=\infty$ (solid green line) and ${\mathcal{N}}_{u}=25$ (solid magenta line) — $\rho_{\otimes}^{-1}$ being proportional to the characteristic length $L$ of the network. In both cases the simulations run for $50000$ generations on a $5000^{2}$ grid. The results were averaged over $25$ simulations. Figure 6 shows that in the absence of death by starvation (${\mathcal{N}}_{u}=\infty$), and for t\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{>}}50, the characteristic length scale of the network grows with time as $\langle L\rangle\propto\langle\rho_{\otimes}^{-1}\rangle\propto t^{1/2}$. On the other hand, in the ${\mathcal{N}}_{u}=25$ case the network attains a regime in which $\langle\rho_{\otimes}\rangle\sim\rm const$, indicating that the average characteristic length scale $\langle L\rangle$ of the network becomes roughly constant in time for t\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{>}}200.

Figure 7 shows the phase space evolution for a single realization of the standard (${\mathcal{N}}_{u}=\infty$) and starvation (${\mathcal{N}}_{u}=25$) models with two and three species (top and bottom, panels, respectively panel). The initial conditions are the same for both standard and starvation models and, in both cases, the simulations run for 10000 generations in a $250^{2}$ grid. Figure 7 confirms that long lasting coexistence only occurs in the starvation model (${\mathcal{N}}_{u}=25$).

Far from the borders of sufficiently large spatial domains occupied by a single species individuals die at each predation attempt, independently of the value of ${\mathcal{N}}_{u}$. On the other hand, an individual is born at each successful reproduction attempt (an unsuccessful predation attempt will occur whenever the passive is an occupied site). Assuming that the probability of the passive being an empty site is equal to $1-\rho_{D}$, where $\rho_{D}$ is equal to the density of individuals deep inside the domains, the equilibrium condition may be written as

[TABLE]

or, equivalently,

[TABLE]

According to Eq. (5) equilibrium is only possible if $p<r$. If $p>r$ then $\rho_{D}=0$.

In practice, the probability of a site far from the borders being occupied is not independent of the distribution of individuals in neighbouring sites, and Eq. (5) is only approximately valid. Fig. 8 shows the comparison between the values of $\rho_{D}$ obtained using the analytical approximation given by Eq. (5) (dashed black line) with the corresponding numerical results (blue dots connected with a solid line) for various values of $p/r$ (assuming that $m=0.5$ and $p+r+m=1$). Fig. 8 shows that the analytical fit is almost perfect for small values of $p/r$ ($p/r<0.5$), but that the numerical results deviate from the analytical prediction for larger values of $p/r$. In particular, the maximum value of $p/r$ consistent with $\rho_{D}\neq 0$ ($p/r\sim 0.75$) is slightly smaller than the analytical prediction ($p/r=1$).

Fig. 9 shows the probability of extinction $P_{\rm ext}$ of at least one species as a function of $p/r$ for ${\mathcal{N}}_{u}=2$ (green dots connected with a solid line) and ${\mathcal{N}}_{u}=25$ (magenta squares connected with a solid line). The results were obtained from an average of $1000$ $200^{2}$ simulations with two species ($N=2$) starting from random initial conditions, again assuming that $m=0.5$ and $p+r+m=1$ (we verified that the probability profiles for $N=3$ are almost identical to those obtained for $N=2$). Fig. 9 shows that the death by starvation, if $p/r$ is sufficiently low, contributes to the preservation of coexistence. The similarity of the probability profiles obtained with ${\mathcal{N}}_{u}=2$ and ${\mathcal{N}}_{u}=25$ (except for a small shift in $p/r$) indicates that the main qualitative results obtained in this paper do not depend on the specific choice of ${\mathcal{N}}_{u}$. For fixed values of $m$, $r$ and $p$, the asymptotic characteristic length scale of the network is larger for larger values of ${\mathcal{N}}_{u}$ (note that for ${\mathcal{N}}_{u}=\infty$, the value of $L$ would always keep growing). Therefore, $L$ becomes larger for ${\mathcal{N}}_{u}=25$ than for ${\mathcal{N}}_{u}=2$. This is the reason why the interval of $p/r$ allowing for coexistence is broader for ${\mathcal{N}}_{u}=2$ than in the ${\mathcal{N}}_{u}=25$ case for simulations performed in a finite box.

5 CONCLUSIONS

In this letter we have investigated the dynamics of spatial stochastic May-Leonard models in two spatial dimensions, considering the death of individuals by starvation after a given number of successive unsuccessful predation attempts as a new ingredient. We have considered models with mutual predation interactions between between any two individuals of different species which, in their standard version, generally lead to the loss of coexistence in a finite time. By considering numerical simulations with two and three species, as well as analytical arguments, we have demonstrated that death by starvation can play an important role on the dynamics of population networks. In particular, we have shown that, if the reproduction rate is sufficiently high, death by starvation may lead be responsible for the preservation of coexistence.

Acknowledgements.

P.P.A. acknowledges the support by FEDER—Fundo Europeu de Desenvolvimento Regional funds through the COMPETE 2020—Operational Programme for Competitiveness and Internationalisation (POCI), and by Portuguese funds through FCT - Fundação para a Ciência e a Tecnologia in the framework of the project POCI-01-0145-FEDER-031938 and the FCT grant UID/FIS/04434/2013. B.F.O thanks Fundação Araucária, Fapern, FCT and INCT-FCx (CNPq/FAPESP) for financial and computational support.