Predator–Prey Dynamics and Ideal Free Distribution in a Heterogeneous Environment

Abstract

The concept of an ideal free distribution (IFD) is extended to a predator–prey system in a heterogeneous environment. We consider reaction–diffusion–advection equations which describe the evolution of spatial distributions of predators and prey under directed migration. Modification of local interaction terms is introduced, if some coefficients depend on resource. Depending on coefficients of local interaction, the different scenarios of predator distribution are possible. We pick out three cases: proportionality to prey (and respectively to resource), indifferent distribution and inversely proportional to the prey. These scenarios apply in the case of nonzero diffusion and taxis under additional conditions on diffusion and migration rates. We examine migration functions for which there are explicit stationary solutions with nonzero densities of both species. To analyze solutions with violation of the IFD conditions, we apply asymptotic expansions and a numerical approach with staggered grids. The results for a two-dimensional domain with no-flux boundary conditions are presented.

1. Introduction

Dynamics of species and their distributions are the most significant issues in mathematical biology and ecology [1]. Currently, cases of heterogeneous habitat and directional movement of species (taxis) described by diffusion–advection–reaction equations are being investigated [2]. Strategies that provide equilibrium states and ensure species survival for extended periods are crucial for population dynamics. The concept of Ideal Free Distribution (IFD), introduced in [3], enables the development of an evolutionarily stable strategy for the interaction of biological communities. Initially, IFD considered only one species, individuals of which have a complete understanding of the habitat and can freely move to any point of it. It was a purely behavioral concept that did not take into account population dynamics. Later, it was expanded to include competing species [4,5] and predator–prey interactions [6,7]. In [5], Ideal Free Distribution with travel costs was used to describe disadvantage for species to move and cessation of movement between different patches. It can be noted that it is promising to analyze the conditions under which IFD is lost and identify the scenarios that are being realized [8,9].

The study of local predator–prey models with IFD has utilized various assumptions about species migration between two patches [10,11]. The works utilized classical Lotka–Volterra equations, modifying predator’s second-kind Holling response and prey’s logistic growth law. In some predator–prey studies, an IFD strategy was associated to game theory and adaptive dynamics [12,13]. It was stated that IFD provides an evolutionarily stable approach. Study [13] explored the impact of resource heterogeneity in the $30\times 30$ cellular automaton model on local predator–prey system dynamics.

The IFD concept was also developed for models described by reaction–diffusion–advection equations [14,15,16,17]. In [14], the authors considered a system of two competing species with different diffusion strategies and established some important results on persistence, extinction, and coexistence of the competing populations with similar growth laws and dispersal strategies. Partially, they discussed the ideal free dispersal strategy when the density of the species matches the carrying capacity of the resources. For one and two competing species, the IFD is the perfect conformity to its resource distribution $p(x,y)$ [12,18]. An integrodifference equation model was analyzed in [19] to find evolutionarily stable dispersal strategies when there is spatial heterogeneity and seasonal variation in a habitat.

The multistability of ideal free distributions was found for the system of two competing species in [20]. There are a number of papers that spread the concept of the IFD on a predator–prey system [6,7]. To enlarge the IFD concept in the case of reaction–diffusion–advection equations in a predator–prey system, it is important to correctly define the growth function of the prey and the functional response of the predator as well as migration phenomena [21].

The goal of this paper is to extend the ideal free distribution concept to a predator–prey system in a heterogeneous environment modeling by reaction–diffusion–advection equations. To reach this aim, we introduce a modified functional response and a growth function for prey. Apart from commonly used predator taxis on prey, we also employ directed migration of prey to the resource and away from predator.

We organize the rest of this paper as follows. The equations and special form of local interaction terms are given in Section 2. Consideration of the system without fluxes is presented in Section 3. Conditions on diffusion and migration parameters supplying ideal free distribution are formulated in Section 4. Numerical simulations and related comments are presented in Section 5. Section 6 deals with the violation of conditions providing ideal free distribution. Theoretical and numerical analysis is presented. A brief summary is given in the last section.

2. Governing Equations for Predator–Prey Dynamics in a Heterogeneous Habitat

A mathematical model of the spatio-temporal interaction of prey with density $u(x,t)$ and predator $v(x,t)$ can be written as a system of equations:

$$\begin{array}{cc}\hfill & \frac{\partial u}{\partial t}=-\nabla q_1+F_1,q_1=-k_1\nabla u+{\alpha }_{1}u\nabla Q_1,-{\alpha }_{2}u\nabla Q_2,\nabla =\left(\frac{\partial }{\partial x},\frac{\partial }{\partial y}\right),\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill & \frac{\partial v}{\partial t}=-\nabla q_2+F_2,q_2=-k_2\nabla v+{\alpha }_{3}v\nabla Q_3.\hfill \end{array}$$

(2)

Fluxes $q_i$, ($i=1,2$) contain the terms corresponding to homogeneous diffusion with coefficients $k_i$ and functions $Q_j$ ($j=1,2,3$), which characterize directed migration with coefficients ${\alpha }_i$.

We note that diffusion terms in the fluxes can be interpreted as “random wandering” of individuals in the area in search of food. However, we mean that diffusion describes a more general process of the natural tendency of each species to spread over the entire range. As for the functions of directed migration, we consider different taxis mechanisms, so that $Q_1=Q_1\left(p\right)$ characterizes the prey resource-oriented taxis (search for the most effective place in terms of food and survival), and $Q_2=Q_2\left(v\right)$ is the prey taxis from danger (the prey try to avoid places where predators congregate). Function $Q_3=Q_3\left(u\right)$ is responsible, respectively, for the taxis of the predator on its resource—the prey.

Functions $F_1$, $F_2$ describe local interactions

$$\begin{array}{cc}\hfill & F_1={\mu }_{1}u{f}_{nm}(u,p)-{\lambda }_{1}g(u,v,p),\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill & F_2={\mu }_{2}g(u,v,p)-{\lambda }_{2}v,\hfill \end{array}$$

(4)

where ${\mu }_i$ and intrinsic growth parameters, ${\lambda }_i$ are mortality coefficients.

In order to expand the description of local interaction and enrich the arsenal of the system’s behavior, we consider the modified prey growth function and the functional response of the predator, introducing a dependence on a generalized resource.

Function $f_{nm}(u,p)$ is defined as

$$f_{n,m}(u,p)=\frac{u^n}{p^m}\left(1-\frac{u}{p}\right),n,m\in \mathbb{N},$$

(5)

and we can describe the logistic growth of the prey ($n=0$, $m=0$), hyperbolic law ($n=1$, $m=0$) and so on. Functional response $g(u,v,p)$ is given in the form of

$$g(u,v,p)=\frac{uv}{p+Cu}.$$

(6)

Parameter $C>0$ allows realization of the Holling functional response of the second kind [22,23]. Modification (6) was proposed in [24] to take into account the non-uniform distribution of the carrying capacity over the environment. The intrinsic growth parameter for prey ${\mu }_1$ and the rate of mortality for predator ${\lambda }_2$ are positive constants as well as coefficients for nonlinear terms ${\lambda }_1$, ${\mu }_2$.

The habitat of both species is a bounded domain, $D=[0,a]\times [0,b]$, with boundary $\partial D$ on which no-flux conditions are set:

$$\begin{array}{c}\hfill q_j·n=0,j=1,2,(x,y)\in \partial D.\end{array}$$

(7)

System (1)–(7) is supplemented with initial distributions of population densities

$$u(x,y,0)=u^0(x,y),v(x,y,0)=v^0(x,y).$$

(8)

3. Local Interaction (Main System without Fluxes)

Let us consider the problem without diffusion and taxis ($k_j=0$, ${\alpha }_j=0$). For brevity, we take $C=0$ and obtain the system of ODEs,

$$\begin{array}{cc}\hfill & \frac{\partial u}{\partial t}=u\left[\frac{{\mu }_1{u}^n}{p^m}\left(1-\frac{u}{p}\right)-\frac{{\lambda }_{1}v}{p}\right]\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill & \frac{\partial v}{\partial t}=v\left[\frac{{\mu }_{2}u}{p}-{\lambda }_2\right].\hfill \end{array}$$

(10)

In this system, spatial coordinates x, y are being parameters because $p=p(x,y)$. There exist several stationary solutions: unstable trivial equilibrium $u=v=0$, solution without predator $u=p$, $v=0$ and coexistence of both species,

$$u=Ap,v=B_n{p}^{\nu },$$

(11)

where

$$A=\frac{{\lambda }_2}{{\mu }_2},B_n=\frac{{\mu }_1}{{\lambda }_1}{A}^n(1-A),\nu =n+1-m.$$

(12)

Solution (11) shows that the prey distribution is proportional to the resource, $p(x,y)$, and does not depend on parameters n and m. The predator distribution is proportional to the degree of resource function $p(x,y)$; order $\nu$ directly depends on n and m.

The following preposition takes place:

Proposition 1.

When $m=n$, Systems (9) and (10) have Stationary solutions (11) and (12) satisfying the concept of ideal free distribution (IFD) for both species.

Proof. 

Indeed, from (11), we obtain $u=Ap$ and $v=Bp$, where constants A and B follow (12). Thus, prey is proportional to resource function $p(x,y)$ and predator is proportional to prey. Therefore, both species correspond to IFD [3,12,14]. □

Corresponding to (11) distribution for any m, n may be named an IFD–like solution because prey is always proportional to resource function $p(x,y)$. Further, we call these solutions IFD for brevity.

Let us discuss the stability of equilibria for Systems (9) and (10). The trivial equilibrium is unstable, solution $u=p$, $v=0$ is stable when ${\lambda }_2\ge {\lambda }_2^{cr}$ and Equilibrium (11) is stable when

$$\left(\frac{n}{n+1}\right){\mu }_2={\lambda }_2^{osc}\le {\lambda }_2<{\lambda }_2^{cr}={\mu }_2.$$

(13)

Here, critical value ${\lambda }_2^{osc}$ corresponds to oscillatory instability, so a limit cycle exists when ${\lambda }_2<{\lambda }_2^{osc}$. The range of stability of stationary solution (11) $\delta \mathsf{\Lambda }$ varies under change in n,

$$\delta \mathsf{\Lambda }={\lambda }_2^{cr}-{\lambda }_2^{osc}=\frac{{\mu }_2}{n+1}.$$

(14)

Different stability conditions can exist for the same $\nu$; as n increases, the stability interval becomes smaller. As a result, we can observe that the prey may control the growth function by particular n and m in order to shorten or lengthen the stability interval for the solution that represents the coexistence of two species. This situation is achieved for the parameter values that fall into the range between the oscillatory regime (${\lambda }_2<{\lambda }_2^{osc}$) and the solution without a predator (${\lambda }_2>{\lambda }_2^{cr}$).

We can obtain the same value of $\nu$ with different combinations of parameters m and n, but parameter n also affects the size of the predator population according to formula

$$B_n=\frac{{\mu }_1}{{\lambda }_1}{\left(\frac{{\lambda }_2}{{\mu }_2}\right)}^n\left(1-\frac{{\lambda }_2}{{\mu }_2}\right),$$

(15)

so that when n increases by N, the ratio of the current value of the predator population to the previous one changes by value

$$\delta B=\frac{B_{n+N}}{B_n}={\left(\frac{{\lambda }_2}{{\mu }_2}\right)}^N,N\in \mathbb{N}.$$

(16)

Since, for a solution with two species, ${\lambda }_2<{\mu }_2$, the predator population decreases. In this sense, the prey, by adjusting n, can not only change the stability interval, but also directly reduce the number of the predator.

4. Ideal Free Distribution and Its Advancement

Consider saving Solutions (11) and (12) in System (1)–(6), when fluxes $q_J$ are nonzero ones.

Proposition 2.

Let directed migration functions $Q_j$ be given as

$$Q_1=\ln p,Q_2=\ln u,Q_3=\ln v,$$

(17)

and the following relations hold:

$${\alpha }_1=k_1+\nu {\alpha }_2,{\alpha }_3=\nu k_2.$$

(18)

Then, (11) is a solution of system (1)–(8).

Proof. 

Firstly, reaction terms are zero because (12) satisfies (9) and (10). Substitution (11) and (12) to fluxes offers

$$\begin{array}{cc}\hfill & q_1=-k_{1}A\nabla p+{\alpha }_{1}Ap\nabla \left(\ln p\right)-{\alpha }_{2}Ap\nabla \left(\ln B{p}^{\nu }\right),\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill & q_2=-k_2\nabla \left(B{p}^{\nu }\right)+{\alpha }_{3}B{p}^{\nu }\nabla \left(Ap\right).\hfill \end{array}$$

(20)

After simplifications, we have

$$\begin{array}{cc}\hfill & q_1=A\nabla p(-k_1+{\alpha }_1+{\alpha }_2{p}^{\nu -1}),\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill & q_2=B{p}^{\nu -1}\nabla p(-k_2\nu +{\alpha }_3).\hfill \end{array}$$

(22)

Using (18), we come to zero fluxes and, respectively, IFD. □

So, with the logarithmic form of directional migration functions (17) and conditions (18), a series of IFD solutions (11) arises, specified by parameter $\nu$.

Table 1. Parameter conditions for IFD solutions.

$\mathit{\nu }$	Conditions
2	${\alpha }_1=k_1+2{\alpha }_2$	${\alpha }_3=2k_2$
1	${\alpha }_1=k_1+{\alpha }_2$	${\alpha }_3=k_2$
0	${\alpha }_1=k_1$	${\alpha }_3=0$ (any $k_2$)
−1	${\alpha }_1=k_1-{\alpha }_2$	${\alpha }_3=-k_2$

presents the conditions on diffusive and migration parameters for several values of parameter $\nu$. It illustrates distinct predator distribution scenarios: proportionality to prey density, neutrality case, and reverse proportionality.

Let us comment on the last two lines, namely cases $\nu =0$ and $\nu =-1$. Biologically, $\nu =0$ corresponds to the absence of predator taxis regarding prey (${\alpha }_3=0$). In the case ${\alpha }_3=-k_2$ ($\nu =-1$), it means that predator taxis is directed to minimum values of prey density. We discuss these scenarios through computational experiments in Section 5.

5. Computer Experiment with IFD and Its Extension

To solve Problem (1)–(8), we apply the method of Lines with staggered grid discretization (see Appendix A) and MATLAB. We begin with parameter values corresponding to $\nu =1$ and the hyperbolic growth law of the prey (weak Allee effect). The calculations use the non-uniform distribution of the prey resource $p(x,y)$, shown in Figure 1. Further, for convenience of analyzing the results, we consider four sections in the areas of local maxima, shown in Figure 1.

Figure 2 shows a three-dimensional picture of the distribution of resource (green), prey (blue) and predator (red). We see that population densities are proportional to the resource function throughout the two-dimensional habitat. The surfaces of the corresponding species lie over each other without intersection. When $\nu =2$, the predator has distribution proportional to the square of the resource; see Figure 3a. This figure shows a comparison of truncated three-dimensional distributions of predator and prey, presented for greater clarity without respect to scale in spatial coordinates. It is clearly visible that the distribution of predator (red surface) is more elongated in places of local maxima compared to that of the prey (blue). We note that if Condition (13) is violated, the system switches into an oscillating mode. Series of Figure 3a–c are shown, where we see how, under the influence of a small disturbance, the system leaves the stationary state. The amplitude increases and the system reaches an oscillatory mode with a constant amplitude; see Figure 3d.

Figure 4a shows the change in the population size of the predator (red color) in comparison with the resource (green color) and the prey population (blue color) for two values of parameter n. Figure 4b illustrates dependencies on parameter n for two characteristics: $\delta B$ (see (16)) and $\delta \mathsf{\Lambda }$ (see (14)) for $N=n$.

Figure 5 shows the distribution of predator and prey at $\nu =0$. Here, there is a constant value of predator population density, independent of the coordinates (red plane), while prey population density is proportional to the distribution of its resource p(x,y) (blue surface). We note that in this case we can also observe an oscillating mode, which occurs when Condition (13) is lost. The same is valid for all rows in Table 1 when $n\ne 0$.

Figure 6 demonstrates a situation where the distribution of predator is inversely proportional to the distribution of prey ($\nu =-1$). The maximum density of the prey corresponds to the minimum for the predator.

6. Violation of Ideal Free Distribution

Let us analyze how the solution changes when the conditions on the parameters are violated. For brevity, we consider the special case of System (1)–(6): $n=m=0$, ${\lambda }_1=1$, ${\mu }_2=1$, $C=0$, ${\alpha }_2=0$:

$$\begin{array}{l}\frac{\partial u}{\partial t}=\nabla (k_1\nabla u-{\alpha }_{1}u\nabla p)+u\left[{\mu }_1\left(1-\frac{u}{p}\right)-\frac{v}{p}\right],\\ \frac{\partial v}{\partial t}=\nabla (k_2\nabla v-{\alpha }_{3}v\nabla u)+v\left[-{\lambda }_2+\frac{u}{p}\right].\end{array}$$

(23)

So, taxis to resource $p\left(x\right)$ acts for the prey, and the predator is oriented only on prey distribution (see case $\nu =1$ in Table 1). We let

$$k_1={\alpha }_1+\epsilon ,k_2={\alpha }_3+\gamma \epsilon .$$

(24)

We look for a stationary solution in the form of a formal series,

$$\begin{array}{ccc}\hfill & & u=u_0+\epsilon u_1+{\epsilon }^2{u}_2+\dots ,\hfill \end{array}$$

(25)

$$\begin{array}{ccc}& & v=v_0+\epsilon v_1+{\epsilon }^2{v}_2+\dots .\hfill \end{array}$$

(26)

When conditions

$$k_1={\alpha }_1,k_2={\alpha }_3$$

are held and $\epsilon =0$, there exists the IFD solution,

$$u_0={\lambda }_{2}p,v_0={\mu }_1(1-{\lambda }_2)p.$$

Collecting terms at the same powers, we obtain the following equations at the first power, $\epsilon$,

$$\begin{array}{ccc}& & 0=\nabla \left[{\alpha }_1\nabla u_1+{\lambda }_2\nabla p-{\alpha }_2\frac{u}{p}\nabla p\right]-{\lambda }_2({\mu }_1{u}_1+v_1)\equiv g_1,\hfill \\ & & 0=\left[{\alpha }_3\nabla v_1+\gamma {\mu }_1(1-{\lambda }_2)\nabla p-{\alpha }_1{\mu }_1(1-{\lambda }_2)p\nabla \left(\frac{u_1}{{\lambda }_{2}p}\right)-{\alpha }_2\frac{v_1}{p}\nabla p\right]+{\mu }_1(1-{\lambda }_2)u_1\equiv g_2.\hfill \end{array}$$

Expression $g_j$ can be written in the form

$$\begin{array}{ccc}\hfill & & g_1=\nabla \left[{\alpha }_{1}p\nabla \left(\frac{u_1}{p}\right)+{\lambda }_2{p}^{\prime }\right]-{\lambda }_2(\mu u_1+v_1),\hfill \end{array}$$

(27)

$$\begin{array}{ccc}& & g_2=\nabla \left[{\alpha }_{3}p\nabla \left(\frac{{\lambda }_1{v}_1-{\mu }_1(1-{\lambda }_2)u_1}{{\lambda }_{2}p}\right)+\gamma {\mu }_1(1-{\lambda }_2)\nabla p\right]+{\mu }_1(1-{\lambda }_2)u_1.\hfill \end{array}$$

(28)

For domain D, from the stationary condition of solution $g_j=0$, integration yields the following integral identities:

$${\int }_D{u}_{1}dx=0,{\int }_D{v}_{1}dx=0.$$

(29)

This means that if the conditions of the IFD are violated, the contribution of first-order additives for deviations from IFD distributions of prey and predator (25) is completely leveled out. This raises the question of whether the additives are in phase, the answer to which is provided by the results of a computational experiment.

Let us analyze these deviations for small variations in diffusion and migration parameters. Figure 7 shows the results of computational experiment at $k_1$ (IFD) and ${\widehat{k}}_1=k_1+\epsilon$ (dashed line). One can see that variations are most noticeable in areas of local maxima of the resource function. With $\nu =0$, the distribution of the predator at $k_1>k_1$ ceases to be homogeneous; we observe the largest deformations in zones of extrema $p\left(x\right)$; see Figure 7c. In the case of $\nu =-1$, Figure 7d, the variations in predator distribution are rather small.

Changing of diffusion coefficient ${\widehat{k}}_2=k_2-\epsilon$ leads to the difference in variations of prey and predator, see Figure 8. For $\nu >0$, there is a decrease in the prey population and an increase in the predator. When $\nu =0$, the picture does not change, as it was promised in Table 1. For negative values of $\nu$, there is an increase in the density of the prey and a decrease in the density of the predator.

7. Discussion and Conclusions

The concept of IFD implies a correlation between the distribution of a species and its resource [3]. For one species, the density is directly proportional to its resource function in a heterogeneous habitat. However, his concept needs to be clarified in the case of a predator–prey system [10]. Our goal is to derive a reaction–diffusion–advection model, being IFD for the prey on a heterogeneous environment, and offer a description of a scenario when a predator occurs on habitat. We can talk about the presence of IFD in the predator–prey system when only the relationship between the prey and its resource is subordinated to this principle, since, in a certain sense, this relationship is primary. Strictly speaking, if the prey is proportional to the resource, which itself is a resource for the predator, then the predator must be proportional to the distribution of the prey to satisfy the concept of IFD. The implementation of stationary solutions in this case is possible for various variants of predator distribution, which are discussed in Section 5.

We design the model under consideration in such a way that both the prey and the predator can influence the dynamics and the stationary distribution of species. For a predator, such a “tool” is the trophic function with parameter C and directed migration function $Q_3$. The prey has a “larger arsenal”: a growth function with parameters m, n and two directed migration functions $Q_1$, $Q_2$. It is stated above that a “pure” IFD solution is possible if the statements of Propositions 1 and 2 are fulfilled.

If simply Statements (17) and (18) of Proposition 2 are met, then only the prey satisfies the IFD concept, and different stationary distributions are possible for the predator. Coefficients m, n greatly impact the shape of stationary distribution as well as the stability interval. We leave open the question of whether it is beneficial for the prey to shorten the stability interval, thereby increasing the risk of occurrence of an oscillatory regime. By controlling parameters m, n in the growth function of the prey, it is possible to reach different distributions, but we are far from thinking that the stationary solutions obtained in this way are the volitional choice of the prey itself. Really, each decision should be viewed as some kind of compromise or superposition of all the acting species. If the predator is immobile—for example, the web of a spider is motionless relative to the free movements of the prey—then the resulting distribution depends only on the prey; see, for example, Figure 5). However, for any behavior of the prey, varying the parameters of diffusion and directed migration, the predator may prevent the scenario chosen by the prey.

The paper attempts to generalize the concept of ideal free distribution to the case of a predator–prey system in a heterogeneous habitat. We extend the concept of ideal free distribution to the case where there is only a proportional distribution of the prey to its resource, while the predator can have different scenarios. We introduce a modernized prey growth Function (5), allowing the implementation of various variants of stationary solutions with an ideal free distribution.

Here, we do not consider the influence of parameter C in (6) on dynamics of the system. This question and temporal dependence may be a subject of future study.