Pattern Formation and Bistability in a Generalist Predator-Prey Model

Abstract

Generalist predators have several food sources and do not depend on one prey species to survive. There has been considerable attention paid by modellers to generalist predator-prey interactions in recent years. Erbach and collaborators in 2013 found a complex dynamics with bistability, limit-cycles and bifurcations in a generalist predator-prey system. In this paper we explore the spatio-temporal dynamics of a reaction-diffusion PDE model for the generalist predator-prey dynamics analyzed by Erbach and colleagues. In particular, we study the Turing and Turing-Hopf pattern formation with special attention to the regime of bistability exhibited by the local model. We derive the conditions for Turing instability and find the region of parameters for which Turing and/or Turing-Hopf instability are possible. By means of numerical simulations, we present the main types of patterns observed for parameters in the Turing domain. In the Turing-Hopf range of the parameters, we observed either stable patterns or homogeneous periodic distributions. Our findings reveal that movement can break the effect of hysteresis observed in the local dynamics, what can have important implication in pest management and species conservation.

1. Introduction

The relevance of considering species movement and spatial distribution is nowadays well known and recognized in mathematical ecology [1]. In some cases, such as the study of biological invasions of an exotic species or the spread of epidemics, the phenomena is essentially spatial. In other situations, however, the explicit introduction of space can change the forecasts of the nonspatial counterpart model. A classical example where the spatial model change the conclusions of the local one, is dynamics of hosts and parasitoids whose coexistence is extensively observed in nature: the nonspatial Nicholson-Bailey host-parasitoid model [2] forecasts extinction of one or both species; however when space and local movement are considered along with this dynamics, the persistence of hosts and parasitoids is obtained and different spatial configurations are observed [3,4].

There are several ways of including space in the model, depending on how the variables are considered, discrete or continuous. If the spatial and temporal variables are assumed continuous and the population is described in terms of densities, Partial Differential Equations (PDE) are the common framework. The so called reaction-diffusion models formulated in terms of second-order PDE have now a long history in theoretical ecology since the analogy between the movement of molecules and the random dispersal of individuals of a population was proposed by Skellam (1951) [5,6,7,8]. Travelling waves of population invasions [9,10,11,12] for studying the spread of epidemics and exotic species; the critical patch size problems which study the minimum region necessary for a population to persist [5,13] and pattern formation in ecology [8,14,15,16] are among the relevant biological phenomena that have been being addressed through reaction-diffusion models.

Biological populations are rarely homogeneously distributed in space what makes the study of pattern formation a very important issue in ecology. Briefly speaking, Turing showed in his seminal paper [17] that, under appropriate conditions, a stationary stable homogeneous distribution of two chemical substances in the absence of diffusion can be destabilized by heterogeneous perturbations when diffusion is present. In their pioneer work Segel and Jakson (1972) [16] first presented the Turing’s ideas of instability induced by diffusion in an ecological context while and Gierer and Meinhardt (1972) [18] identified the “activator-inhibitor mechanism”. Since then the literature abounds with research on the Turing mechanism of pattern formation in Biology [19] and, in particular, in predator-prey systems. When the stable equilibrium of the local dynamics undergoes a Hopf bifurcation as some parameter varies, a limit cycle, very common in predator-prey systems, appears. For parameters close to the bifurcation point, it is known that heterogeneous perturbations can lead the spatial system to complex spatio-temporal dynamics including spatio-temporal chaos. This mechanism of pattern formation is known as Turing-Hopf bifurcation [20,21].

Models focusing on pattern formation in predator-prey systems, in most cases, consider that the predator is a specialist [20,21,22]. It means that predators feed itself on a single species of prey, so that in the absence of it the predator goes extinct. Generalist predators, in their turn, have more than one food source and do not depend exclusively on one prey to survive. Although generalist predators have been receiving considerable attention in the context of invasions and biological control in the last decades, there are few studies analysing the influence of a generalist predator in a spatially distributed system [23,24]. A recent study analyse how a generalist predator may influence the stable heterogeneous spatial pattern formation [25]. In this research, the author considered a reaction-diffusion generalist predator-prey system with constant alternative food source and Holling type II functional response [26]. Chakraborty (2015) [25] analysed the case in which the model presents just one coexistence equilibrium. As a result of Turing instability, patterns such as cold spots, hot spots and stripes were observed. For parameters near the Turing-Hopf region, the model exhibits a chaotic behavior in space and time.

The model by Erbach and collaborators (2013) [27] consists of a prey species, which grows logistically, and a predator that has a constant alternative food source, so that it does not become extinct in the absence of the focal prey. Predation is described by the functional response Holling type III. The local model explored by Erbach and collaborators shows a rich and complex dynamics which can exhibit one stable interior equilibrium, bistability of two stable equilibria, limit cycles (stable and unstable) and bistability between a stable equilibrium and a stable limit-cycle with plenty of bifurcations. The concomitant existence of two stable equilibria creates basins of attraction which means that the initial condition determines the fate of the system solution. It also makes possible sudden and abrupt state transitions that may be dramatic from the ecological, social or economic point of view [28].

Erbach and collaborators (2013) say, in the discussion of the above mentioned paper, that they were taken by a great curiosity about the effects of diffusion on the spatiotemporal dynamics of the species conducted by this dynamics. Hence, here we explore the effects of spatial distribution and movement on the generalist predator-prey dynamics studied by Erbach et al. (2013). Our goal is to investigate the role of space on the complex dynamics of their model. In particular, we intend to analyze the Turing and Turing-Hopf pattern formation with special attention to the effects of dispersal on the bistability regime. We start describing the model and giving a brief review of the main results by Erbach and collaborators in Section 2. In what follows, we derive the Turing instability conditions for diffusive instability in the model in Section 3 and then, by means of numerical simulations, presented in Section 4, we illustrate the patterns obtained. Finally, we discuss the implications of our results for the dynamics of the species.

2. The Model

The proposed model is a reaction-diffusion system with a prey and a generalist predator in the presence of a constant alternative food for the predator. Predation is described by the functional response Holling type III, which means that the rate in which predators attack the prey is small when the prey density is low. In this scenario, the predator prefers to capture more abundant alternative prey species and abandon the focal prey. However, as the focal prey density increases, the rate at which predators capture it, also increases, in a phenomenon known as prey switching. Since for many biological species the individuals disperse randomly in space [8], we adopt the classical diffusion as paradigma for modelling the spatial dispersal. Hence, the spatio-temporal dynamics is described by the following system:

$$\begin{array}{c}{\displaystyle \frac{\partial N}{\partial T}}=RN\left(1-{\displaystyle \frac{N}{K}}\right)-{\displaystyle \frac{A{N}^{2}P}{1+H{N}^2}}+D_1\left({\displaystyle \frac{{\partial }^{2}N}{\partial X^2}}+{\displaystyle \frac{{\partial }^{2}N}{\partial Y^2}}\right),\hfill \\ {\displaystyle \frac{\partial P}{\partial T}}={\displaystyle \frac{BP}{1+EP}}-MP+C{\displaystyle \frac{A{N}^{2}P}{1+H{N}^2}}+D_2\left({\displaystyle \frac{{\partial }^{2}P}{\partial X^2}}+{\displaystyle \frac{{\partial }^{2}P}{\partial Y^2}}\right),\hfill \end{array}$$

(1)

where $N(X,Y,T)$ and $P(X,Y,T)$ denote the densities of the prey and the predator, respectively, at $(X,Y)\in {\mathbb{R}}^2$ and time $T\ge 0$.

In the absence of predator, the prey population reproduces according to the logistic function. Parameters R and K denote the intrinsic growth rate and the carrying capacity of the prey. In its turn, in the absence of the focal prey, the predator population reproduces following a Beverton-Holt function. The parameter B denotes the per capita reproduction rate of the predator while E represents the strength of density-dependence. The death rate of the predator is denoted by M. In the term corresponding to the type-III functional response, A represents the searching efficiency and H, the handling time [26]. The efficiency of converting prey into predator biomass is denoted by C. $D_1$ and $D_2$ are the diffusion coefficients of prey and predator, respectively. All the parameters are positive and since we are considering that predator can survive in the absence of the focal prey, we assume $B>M$.

We introduce the following dimensionless variables: $n=N/K$, $p=\left(AKP\right)/R$, $x=X\sqrt{R/D_2}$, $y=Y\sqrt{R/D_2}$ and $t=RT$, so that predator and prey population dynamics are governed by the nondimensional system

$$\begin{array}{c}{\displaystyle \frac{\partial n}{\partial t}}=n(1-n)-{\displaystyle \frac{n^{2}p}{1+a{n}^2}}+D\left({\displaystyle \frac{{\partial }^{2}n}{\partial x^2}}+{\displaystyle \frac{{\partial }^{2}n}{\partial y^2}}\right),\hfill \\ {\displaystyle \frac{\partial p}{\partial t}}={\displaystyle \frac{cp}{1+dp}}-ep+{\displaystyle \frac{b{n}^{2}p}{1+a{n}^2}}+\left({\displaystyle \frac{{\partial }^{2}p}{\partial x^2}}+{\displaystyle \frac{{\partial }^{2}p}{\partial y^2}}\right),\hfill \end{array}$$

(2)

where $a=H{K}^2$, $b=\left(CA{K}^2\right)/R$, $c=B/R$, $d=\left(ER\right)/\left(AK\right)$, $e=M/R$ and $D=D_1/D_2$. Model (2) now contains only six parameters (against ten in its dimensional form). Since $B>M$, we have now $c>e$.

Let $\mathsf{\Omega }$ be a two-dimensional bounded square domain and $\eta$ the outward unit normal vector of the boundary $\partial \mathsf{\Omega }$. We investigate the model (2) under the following initial conditions:

$$\begin{array}{c}{n}_0=n(x,y,0)>0,\hfill \\ p_0=p(x,y,0)>0,\hfill \end{array}$$

(3)

where $(x,y)\in \mathsf{\Omega }$, with zero-flux boundary conditions

$${\displaystyle \frac{\partial n}{\partial \eta }}{|}_{(x,y)}={\displaystyle \frac{\partial p}{\partial \eta }}{|}_{(x,y)}=0,$$

(4)

where $(x,y)\in \partial \mathsf{\Omega }$.

Local Stability

For completeness, we briefly summarize the main results of the local stability. A detailed analysis can be found in [27].

The corresponding non-diffusive model is

$$\begin{array}{c}{\displaystyle \frac{dn}{dt}}=n(1-n)-{\displaystyle \frac{n^{2}p}{1+a{n}^2}},\hfill \\ {\displaystyle \frac{dp}{dt}}={\displaystyle \frac{cp}{1+dp}}-ep+{\displaystyle \frac{b{n}^{2}p}{1+a{n}^2}}.\hfill \end{array}$$

(5)

Moreover $n=0$ and $p=0$, the n-nullclines and p-nullclines are given, respectively, by

$$f\left(n\right)={\displaystyle \frac{(1+a{n}^2)(1-n)}{n}}\mathrm{and}g\left(n\right)={\displaystyle \frac{1}{d}}\left({\displaystyle \frac{c(1+a{n}^2)}{e(1+a{n}^2)-b{n}^2}}-1\right).$$

(6)

The shape of function f depends on parameter a. If $a\le 27$ then $f\left(n\right)$ is monotonically decreasing, whereas if $a>27$, $f\left(n\right)$ has two local extrema ([27]; see also [22]). Moreover, the n-nullcline always has a singularity at $n=0$ and a root at $n=1$. Function $g\left(n\right)$, on the other hand, is always monotonically increasing. When $b<ae$, $g\left(n\right)$ is bounded. However, for $b>ae$, there is a vertical asymptote at $n=\sqrt{e/(b-ae)}$ [27]. Figure 1 shows a typical shape of both nullclines for $a>27$ and $b<ae$.

The model always has the trivial equilibrium $(0,0)$ and the two semi-trivial equilibria $(1,0)$ and $(0,(c-e)/ed)$. Furthermore, depending on the parameter values, the system can have one or three coexistence equilibria.

Linearization of the system around the coexistence equilibrium $(n^{*},p^{*})$ produces the Jacobian matrix

$$J(n^{*},p^{*})=\left(\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right),$$

(7)

where $a_{11}=1-2n^{*}-{\displaystyle \frac{2n^{*}{p}^{*}}{{(1+a{n}^{*})}^2}}$, $a_{12}=-{\displaystyle \frac{{n^{*}}^2}{1+a{n^{*}}^2}}$, $a_{21}={\displaystyle \frac{2b{n}^{*}{p}^{*}}{{(1+a{n^{*}}^2)}^2}}$ and $a_{22}={\displaystyle \frac{b{n^{*}}^2}{1+a{n^{*}}^2}}+{\displaystyle \frac{c}{{(1+d{p}^{*})}^2}}-e$.

By applying the Routh-Hurwitz criterion, the equilibrium $(0,0)$ is always an unstable node and $(1,0)$ and $(0,(c-e)/ed)$ are saddles. To analyze the stability of the coexistence equilibria, Erbach and colleagues [27] found two relations between the slope of f and g with the elements of matrix J, as follows:

$$\begin{array}{ccc}\hfill f^{\prime }\left(n^{*}\right)& =& -a_{11}/a_{12},\hfill \end{array}$$

(8)

$$\begin{array}{ccc}\hfill g^{\prime }\left(n^{*}\right)& =& -a_{21}/a_{22}.\hfill \end{array}$$

(9)

These conditions lead to the following conclusions:

(i)

If $f^{\prime }\left(n^{*}\right)<0$, then the equilibrium is stable;

(ii)

The equilibrium is a saddle point if and only if $f^{\prime }\left(n^{*}\right)>g^{\prime }\left(n^{*}\right)$.

In other words, conclusion (i) establishes that if an equilibrium point is located where the n-nullcline is decreasing, then it is stable. From this, Erbach and colleagues conclude that for $a<27$ the unique coexistence equilibria is always stable. Conclusion (ii) states that if an equilibrium point stands where the slope of the n-nullcline is greater than the slope of the p-nullcline, then this equilibrium is a saddle point. This scenario only occurs for $a>27$ and a proper choice for the other parameters. The authors found a complex dynamical behaviour including biestability, limit cycles and several bifurcations related to the positive slope of n-nullcline at the equilibrium. Pattern formation in the spatial model will be explored in these different scenarios.

3. Turing Instability

The Turing mechanism [17] shows that, under appropriate conditions, two chemical substances that react and diffuse can lead to stable heterogeneous spatial patterns. Turing considered Fickian classical diffusion for the substances movement and the law of mass action for the interaction. For diffusive instability to occur it is assumed that the reaction achieved a stable steady state. The classical diffusion is a macroscopic description (or on the population level) of the random movement of the molecules (or biological individuals) [29]. This process results in a net movement from regions of high concentrations towards regions of low concentrations; it means that it has a stabilizing effect that lead to an uniform spatial distribution. Counter-intuitively, Turing showed that the combination of two stabilizing processes (diffusion and stable interaction) can have a destabilizing effect that can promote heterogeneous distributions of the substances.

In this section, we present the conditions that lead to diffusive instability for model (2), following the usual analysis (see, e.g., [6,15,16,30]). Briefly speaking, Turing instability consists in two assumptions: (1) initially, in the absence of diffusion, the system must exhibit a stationary homogeneous state and (2) in the presence of diffusion, the stationary homogeneous state is unstable to small perturbations, what can lead to spatially stable heterogeneous patterns.

Therefore, let us consider the existence of a stationary uniform solution $E^{*}=(n^{*},p^{*})$ for the system (2). To examine the stability of $E^{*}$ with respect to perturbations

$$\begin{array}{ccc}\hfill \epsilon (x,y,t)& =& {\alpha }_{1}cos(k_{1}x+k_{2}y)e^{\lambda t},\hfill \end{array}$$

(10)

$$\begin{array}{ccc}\hfill \delta (x,y,t)& =& {\alpha }_{2}cos(k_{1}x+k_{2}y)e^{\lambda t},\hfill \end{array}$$

(11)

where ${\alpha }_i$, $k_i$ ($i=1,2$) and $\lambda$ are constants, we linearize system (2) and analyze the eigenvalues $\lambda$ of the matrix

$$J_D=\left(\begin{array}{cc}{a}_{11}-D{k}^2& a_{12}\\ a_{21}& a_{22}-k^2\end{array}\right),$$

(12)

where $k^2=k_1^2+k_2^2$ is the wave number (see [16] for a detailed analysis).

Perturbations (10) and (11) decay if and only if $Re\left(\lambda \right)<0$. In other words, the homogeneous equilibrium $(n^{*},p^{*})$ of system (2) is stable if $tr{J}_D(n^{*},p^{*})<0$ and $det{J}_D(n^{*},p^{*})>0$, i.e.,

$$\begin{array}{ccc}\hfill a_{11}+a_{22}-(D+1)k^2& <& 0,\hfill \end{array}$$

(13)

$$\begin{array}{ccc}\hfill (a_{11}-D{k}^2)(a_{22}-k^2)-a_{12}{a}_{21}& >& 0.\hfill \end{array}$$

(14)

Note that for $k=0$ (which means perturbations of zero wavenumber), inequalities (13) and (14) become

$$\begin{array}{ccc}\hfill a_{11}+a_{22}& <& 0,\hfill \end{array}$$

(15)

$$\begin{array}{ccc}\hfill a_{11}{a}_{22}-a_{12}{a}_{21}& >& 0.\hfill \end{array}$$

(16)

These conditions have already been satisfied since we assume the existence of a stationary homogeneous state in the absence of diffusion.

Now, by violating (13) or (14), the stationary homogeneous state becomes unstable to small perturbations in the presence of diffusion, satisfying the second assumption of diffusive instability. But inequality (13) is always true, because $a_{11}+a_{22}<0$. So inequality (14) is the only one that can be “broken”.

Expanding the left hand side of (14) and reversing the inequality, we find $H\left(k^2\right)$ given by

$$H\left(k^2\right)=D{k}^4-(D{a}_{22}+a_{11})k^2+a_{11}{a}_{22}-a_{12}{a}_{21}<0.$$

(17)

Function $H\left(k^2\right)$ is a quadratic function and its graph is a parabola. Since we are looking for a positive $k^2$ such that $H\left(k^2\right)<0$, it is necessary that

$$D{a}_{22}+a_{11}>0.$$

(18)

Since $D=D_1/D_2$, from (18) follows that $D\ne 1$ (or $D_1\ne D_2$). Furthermore, the minimum of H is given by

$$H\left(k_{min}^2\right)=a_{11}{a}_{22}-a_{12}{a}_{21}-{\displaystyle \frac{{(D{a}_{22}+a_{11})}^2}{4D}},$$

(19)

where $k_{min}^2={\displaystyle \frac{D{a}_{22}+a_{11}}{2D}}$.

In order to satisfy inequality (17) is sufficient that H be negative at its minimum. Therefore, imposing $H\left(k_{min}^2\right)<0$ we obtain the last condition for diffusive instability:

$${(D{a}_{22}+a_{11})}^2-4D(a_{11}{a}_{22}-a_{12}{a}_{21})>0.$$

(20)

In summary, from (15), (16) and (20), the conditions for Turing instability are:

$$\begin{array}{ccc}\hfill a_{11}+a_{22}& <& 0,\hfill \end{array}$$

(21)

$$\begin{array}{ccc}\hfill a_{11}{a}_{22}-a_{12}{a}_{21}& >& 0,\hfill \end{array}$$

(22)

$$\begin{array}{ccc}\hfill {(D{a}_{22}+a_{11})}^2-4D(a_{11}{a}_{22}-a_{12}{a}_{21})& >& 0.\hfill \end{array}$$

(23)

As mentioned above, the Turing instability starts with a homogeneous locally stable steady state in the absence of diffusion what means that the system returns to the spatially homogeneous equilibrium distribution after being homogeneously perturbed. However, if heterogeneous small perturbations of a certain wave number are applied, the system can evolve to heterogeneous spatial distributions. If the homogeneous equilibrium stability is broken via a Hopf bifurcation, the local dynamics becomes oscillatory and homogeneous small perturbations do not vanish with time. The Turing-Hopf bifurcation is then characterized as the point in the parameter space for which the system steady state is destabilized by both homogeneous (Hopf) as well as heterogeneous (Turing) perturbations [20,21,31,32]. Many authors have reported that in the vicinity or inside the Turing-Hopf domain heterogeneous stationary patterns [21] as well as spatially irregular non-stationary patterns and spatiotemporal chaos [20,31] are expected.

4. Numerical Simulations

Our goal in this section is to investigate the effects of species movement on the dynamics, in particular, we intend to investigate the existence of spatial patterns in the system (2). We perform extensive numerical simulations for the model (2) in a bounded spatial habitat $\mathsf{\Omega }=[0,12.5]\times [0,12.5]\subset \mathbb{R}\times \mathbb{R}$ with reflective zero-flux boundary conditions. To solve the system, we discretized the model in time and space using FTCS (forward-time central-space) scheme in a personal code developed in Mathematica 10.0, in which we take $\Delta x=\Delta y=0.25$ and the time-step is $\Delta t=0.01$. From an initial random perturbations of the uniform distribution of the species at the local equilibrium, we run the simulations until the system reach a steady state. In all the cases, predator and prey species assumed the same type of spatial distribution so that we will only refer to the prey one.

The simplest case occurs when we take $a\le 27$ so that f (n-nullcline) is always decreasing and thus the model has only one coexistence equilibrium. In this case, although conditions (21) and (22) are satisfied, inequality (23) never happens. Indeed, since $a\le 27$, the slope of f is negative and from conclusion (i), the equilibrium state is stable. Now, noting that $a_{12}$ is always negative, from Equation (8) we conclude that $a_{11}$ must be negative too. Moreover, as g is always increasing ($g^{\prime }>0$) and $a_{21}>0$, Equation (9) implies that $a_{22}$ is negative. Therefore, as we have $a_{11}<0$, $a_{22}<0$ and $D>0$, condition (23) is not satisfied, which means that there is no diffusive instability for $a\le 27$.

From now on, we focus the analysis on the case $a>27$ and study the following three sets of parameters a, c, d and e for which the dynamics has been explored by Erbach and collaborators (2013):

• Set I: $a=100$, $c=1$, $d=0.9$ and $e=0.5$;
• Set II: $a=100$, $c=1$, $d=0.1$ and $e=0.5$;
• Set III: $a=80$, $c=0.7$, $d=0.6$ and $e=0.1$.

For each set, we present the corresponding bifurcation diagram with b as controlling parameter. After applying the Turing conditions for diffusive instability, we show in the space of parameters $b,D$, the regions for which Turing and/or Turing-Hopf instability can occurs, referred to as Turing or Turing-Hopf region, accordingly.

We also present the density plots of the species distribution to illustrate the patterns obtained. We will refer to the heterogeneous distributions obtained as Turing or Turing-Hopf patterns, accordingly. We adopted the classification used by Baurmann and collaborators [20] for the patterns obtained: Homogeneous Distribution and Stationary Patterns. The first one obviously refers to a homogeneous distribution of the species over the whole habitat. If the system evolves to a stationary heterogeneous spatial distribution of prey and predators, we will use the second classification. The resulting pattern can have different spatial configurations. Isolated areas with high population densities surrounded by areas with very low population densities are called hot spots. However, we make a distinction when high picks of density are observed and call the pattern spikes. On the other hand, when the distribution of the species presents isolated areas with low densities, the resulting pattern is called cold spots. When high and low population densities are alternated, the pattern is called stripes. Patterns that do not fit any of these characteristics will be simply called heterogeneous patterns.

5. Results

5.1. Set I

For parameters’ set I the system presents only one coexistence equilibrium which can lose stability and gives place to a stable limit cycle depending on the values of b. Figure 2a illustrates a bifurcation diagram of the prey equilibrium with respect to b while Figure 2b shows one nullcline for the prey and several predator-nullclines. In both figures the values of b for which we obtain Turing as well as Turing-Hopf spatial patterns can be observed.

From our previous comments, Turing instability can only be possible when the equilibrium occurs in an increasing region of the n-nullcline. For parameters’ set I, this happens for b in the gray interval of Figure 2a. For b in the dark gray region, the equilibrium is stable and Turing instability condition (23) can be satisfied for some values of D. For parameters in light gray interval, the coexistence equilibrium is unstable and limit cycles are observed what suggests possible Turing-Hopf instabilities for suitable values of D.

Through this analysis we restricted our attention to the range of b for which diffusion driven instabilities are possible. We then, numerically verify the values of D for which condition (23) is satisfied. The only restriction for D is $0<D<1$ since $D_2<D_1$. Figure 3 illustrates the region in the $b,D$ parameters space where Turing and/or Turing-Hopf instabilities take place while Figure 4 indicates the types of pattern obtained for different combinations of the parameters b and D inside these regions.

Examples of patterns obtained for parameters inside the Turing-Hopf region are illustrated in Figure 5. We can observe that increasing b and keeping $D=0.01$, the resulting pattern evolves from cold spots to hot spots after passing by a transition of stripes patterns. For grater values of D even though inside the Turing-Hopf region, the system evolves to a homogeneous population distribution. Since the equilibrium of the local dynamics is not stable for the parameters Set I with b inside this region; the system undergoes oscillations according to the corresponding limit cycle. Hence, the spatial dynamics exhibit time-oscillating homogeneous distributions.

For parameters taken inside the Turing instability region of the parameters illustrated in Figure 4b we obtained cold spots as resulting patterns (see Figure 6a,b for examples). On the other hand, with parameters inside the Turing region illustrated in Figure 4c, spikes and hot spots were observed as, for example, those showed in Figure 6c,d, respectively.

5.2. Set II

For parameters set II, the system exhibits bistability and unstable limit cycles as illustrated in Figure 7. It suggests that it would be possible, for the same parameter range, two different pattern regimes depending on whether the initial perturbation is made around one or another equilibrium. Furthermore, one can ask about the influence of one equilibrium on the other. We will refer to the equilibrium with small prey density as $E_1$ and to the equilibrium with high prey density as $E_2$.

For $b<20.23$, $E_2$ is located where the n-nullcline is decreasing. Hence, condition (23) is not satisfied and Turing pattern is not possible for any perturbation of this equilibrium. For $20.92<b<25.03$ approximately, which includes regions $R_2$, $R_3$ and $R_4$ of Figure 7, the diffusive instability conditions (21)–(23) are satisfied for equilibrium $E_1$. For b in the small interval $R_3$ Turing instability is possible for both equilibria.

We then applied condition (23) for both equilibria in order to find the combinations of b and D for which Turing instability can happen. Figure 8a,b illustrate these regions for the equilibrium $E_2$ and $E_1$, respectively. For b in $R_3$ and D taken in the dark gray region of Figure 8a, small perturbation of equilibrium $E_2$, lead the species distributions to cold spots. For the equilibrium $E_1$, on the other hand, we obtained different results. For very small values of D ($D<0.005$ approximately) regardless the value of b in the interval $20.92<b<25.03$ which includes regions $R2$ to $R4$, small heterogeneous initial perturbations lead to non-homogeneous patterns as those illustrated in Figure 9. On the other hand, as D increases, cold spots and heterogeneous stripes like patterns as illustrated in Figure 10 are obtained in regions $R2$ and $R3$.

An interesting phenomenon is observed in region $R2$. For b close to region $R_3$, initial perturbations of equilibrium $E_1$ lead to heterogeneous cold spots as those illustrated in Figure 5a. However, for small values of b in region $R_2$ (approximately $20.92<b<22$) and $D\ge 0.005$ initial perturbations around $E_1$ grow and the prey is lead to a homogeneous distribution at the equilibrium $E_2$ as represented in Figure 11. The blue points indicate the prey and predator densities in all the domain. For $t=0$ (Figure 11a), all the points are around the equilibrium $E_1$. Without diffusion, all the values would converge to $E_1$. However, with diffusion, the solution in all the domain migrate during the transients (see Figure 11b–e) and finally tend to equilibrium $E_2$ (Figure 11f), even though the parameters belong to the Turing region. Nevertheless, there is no contradiction with Turing’s theory since it does not state that when diffusion is present, heterogeneous patterns will necessarily takes place. It says that the homogeneous distribution is unstable to small perturbations, what is the case here: perturbations around $E_1$ grow and the system is led to a uniform distribution at the equilibrium value $E_2$. This phenomenon promotes a rupture of the effect of hysteresis observed in the local dynamics; in the spatially distributed model, the system state jumps from $E_1$ to $E_2$ due to perturbations caused by movement even though the initial distribution is inside the basin of attraction of $E_1$. We also observe that the unstable limit-cycles observed for small values b in $R_2$ (see Figure 7) do not have any influence on the patterns obtained (for $D<0.005$) and on the solution migration to $E_2$ (for greater values of D).

5.3. Set III

Parameters’ set III presents a much more complex local dynamics than the previous two cases. Differently from Set II, no bistability is observed. However, there is the coexistence of a stable equilibrium with a stable limit-cycle. We detected six intervals of b with different behaviour related to pattern formation which will be explained in what follows.

Figure 12 illustrates the bifurcation diagram for the prey with parameter b as controlling parameter. We will call $E_1$ the prey equilibrium with lower density and $E_2$ the higher prey equilibrium value. Figure 13 and Figure 14 show the combinations of parameters b and D for which diffusive instability is obtained, respectively, for equilibrium $E_1$ and $E_2$ in different intervals for b. The dots in these two figures indicate the combinations tested and the type of pattern obtained.

In the intervals corresponding to $R_1$ of Figure 12, $E_1$ is the only stable equilibrium. All the diffusive instability conditions are satisfied and hence, Turing patterns are possible for $E_1$ in $R_1$. For D in gray region of Figure 13c, different heterogeneous patterns and stripes were obtained.

In $R_2$ of Figure 12, $E_1$ is still stable but now it is surrounded by two limit-cycles: one unstable and the other stable. Turing patterns are possible for $E_1$ in $R_2$ since all the conditions for diffusive instability can be satisfied for suitable combinations of b and D. The resulting patterns for D taken inside the corresponding gray region $R_2$ of Figure 13c are heterogeneous patterns. Since for b in $R_2$ there is an unstable limit-cycle surrounding $E_1$, we initially perturbed this equilibrium up to $1\%$ of its value in each point.

For b in region $R_3$, $E_1$ becomes unstable but the stable limit cycle is still surrounding it (see Figure 12). The condition (23) can still be satisfied what characterize region $R_3$ as a region where Turing-Hopf patterns are possible. For small values of D inside the gray region of Figure 13b, heterogeneous patterns and cold spots are obtained. However, with larger values of D, still inside the Turing region, the system is driven to a homogeneous periodic distribution which assumes the corresponding values of the stable limit cycle that surrounds $E_1$.

Regions $R_4$ and $R_5$ are Turing-Hopf regions for $E_1$ since it is an unstable equilibrium surrounded by a stable limit-cycle. The combinations of b and D for which diffusive instability is observed can be seen in Figure 13a. Small perturbations of the homogeneous equilibrium distributions lead to heterogeneous as well as homogeneous periodic solutions for b and D in these regions. Cold spots patterns appear for low values of D in both regions $R_4$ and $R_5$ while homogeneous periodic solutions can be seen for larger values of D.

A saddle-node bifurcation takes place as b cross the border of region $R_3$ towards $R_4$ and the equilibrium $E_2$, at high prey levels, appears. $E_2$, which is also unstable for all values of b in $R_4$, is surrounded by a stable limit-cycle. However, condition (23) holds for b and D in a certain range of values as illustrated in Figure 14c, which characterizes it as a Turing-Hopf region for equilibrium $E_2$. We tested several combinations of parameters in this region and observed cold spots patterns for low values of D. For higher values of D, still inside the gray region, homogeneous periodic solutions are obtained. We can observe an unstable limit-cycle around $E_2$; with high values of D, the population values in all the domain jump outside this cycle and the system is lead to the stable outer cycle. It is worth noting that $R_4$ is also a Turing-Hopf region for equilibrium $E_1$ as mentioned above.

Moving b to region $R_5$, equilibrium $E_2$ becomes stable. However, it is surrounded by an unstable limit-cycle and by the outer stable limit-cycle. All the conditions (21) to (23) can be satisfied for $E_2$ in this region. For b and D the region illustrated in Figure 14b, the Turing patterns obtained are cold spots.

Finally, for b in region $R_6$, both limit-cycles disappear and $E_2$ is still stable. In $R_6$ the diffusive instability conditions can be satisfied only for $E_2$ if D is taken in the appropriate range illustrated in Figure 14a. Cold spots are the kind of distribution we obtained in this case.

6. Discussion

In this paper, we have analysed the effects of spatial distribution and species movement on the generalist predator-prey dynamics studied by Erbach et al. (2013) [27]. In particular, we have investigated the scenarios of pattern formation in the rich dynamics that the above mentioned authors proposed. Our main goal was to elucidate how is pattern formation influenced by limit-cycles and bistability present in the rich local dynamics. We then derived the conditions for Turing instability and developed numerical simulations for different sets of relevant parameters for the local kinetics. We found the regions in the parameter space b (related to the predation intensity) and D (the relation between prey and predator motility) for which Turing and Turing-Hopf patterns can be obtained, in most cases, for parameters far way from the bifurcation. By means of numerical simulations, we observed that hot spots, cold spots, stripes and a mixture of them can emerge depending on the predation intensity. With small values of b, predation on the focal prey is low and we observe isolated regions of lower prey (and predator) density; that is, cold spots are the typical patterns for low values of b. On the other hand, when predation is high (great values of b), the prey remains confined to limited regions of habitat in the so called hot spots patterns if predator movement rate is great compared to the prey. As consequence, the type of pattern can be associated with the effectiveness of the predator on controlling the focal prey. If hot spots are observed, prey density is around its low equilibrium level.

6.1. Spatiotemporal Chaos

In the Turing-Hopf instability domain, we obtained either stationary patterns (for small values of D, including hot/cold spots and stripes depending on b) or oscillating homogeneous distributions (for larger values of D). Spatiotemporal chaos, frequently observed in specialist [20,21] as well as generalist [25] predator-prey systems for parameters, close or inside, the Turing-Hopf instability regions, was not observed for any combinations of parameters, neither close the bifurcation nor into the region where chaotic behaviour would be expected. We also simulated Equation (2), for the parameters set III with $b=4.2$ and $D=0.05$; that is, inside the Turing-Hopf domain, from the initial perturbations $n(x,y,0)=n^{*}+0.01sin\left(\epsilon (y-20)\right)$ and $p(x,y,0)=p^{*}+0.01sin\left(\epsilon (y-20)\right)$, for $\epsilon =0.006$, $\epsilon =0.025$ and $\epsilon =0.25$. In all these simulations the system also converged to either stationary heterogeneous patterns or homogeneous oscillating distributions. None of the initial perturbation tested, produced spatiotemporal chaos.

Contrary to [20,31] who found spatially irregular non-stationary patterns and spatiotemporal chaos for parameters in the vicinity of a Turing-Hopf bifurcation, Banerjee and Petrovskii [21] obtained heterogeneous stationary patterns for parameters close to this type of bifurcation. Spatiotemporal chaos, nonetheless, was obtained in [21] for parameters into the Turing-Hopf range. Furthermore, there are evidences that Turing instability is not necessary to trigger chaos; Banerjee and Petrovskii [21] showed that chaotic spatial patterns can appear in a specialist predator-prey system, outside the Turing domain, as a result of non-homogeneous perturbations in a region of parameters where limit-cycles are present in the local dynamics. We then tested the above mentioned parameters in set III with $D=1$ (that is, outside the Turing domain) and again no irregular spatiotemporal pattern was observed. Chakraborty [25] have obtained chaotic spatiotemporal oscillations in a generalist predator-prey interaction, however in his model, predation occurs according to the Holling type II functional response. Our results, on the other hand, indicate that neither Turing-Hopf nor Hopf bifurcations can always be associated with spatiotemporal chaos in continuous reaction-diffusion predator-prey systems. A similar behaviour has already been found by Rodrigues and colleagues [32] in a discrete Coupled Map Lattice for predator-prey system, in which stationary patterns are observed even when the local dynamics is oscillatory. Although more investigation is needed, our findings suggest that the Holling type III function response such as we have used in the present work, may have a stabilizing effect on the spatiotemporal behaviour of the generalist predator-prey dynamics.

6.2. Bistability and Regime Shifts

The model studied here presents bistability that emerge from saddle-node bifurcations, which is the scenario for jumps and hysteresis to occur [33].

When the prey is an agricultural pest controlled by the generalist predator, the hysteresis effect observed in the local dynamics means that once the prey density falls from the high equilibrium density $E_2$ to the desired low equilibrium $E_1$, as a result of increasing b, predation efforts can be relaxed to keep the control. That is, the pest density would only grow back to the high level $E_2$ if the predation decreases below the leftmost saddle-node bifurcation point, called now $b_{c_1}$ (see Figure 7). Our results, nonetheless, show that spatial movement of the species (as it occurs with many species in real world plantations) may cause a shift back of the prey density to the upper levels if predation intensity is slightly relaxed below the level that causes the prey density decay, the rightmost saddle-node bifurcation point $b_{c_2}$. Close to this parameter value, the basin of attraction of the lower density equilibrium $E_1$ is considerable large (see Figure 7) so that, small perturbations caused by movement do not promote a shift back to high density $E_2$. However, as the parameter b decreases from $b_{c_2}$, the basin of attraction of $E_1$ and consequently, its resilience, diminishes. Now, in the reaction-diffusion model, for values of $b>b_{c1}$, perturbations due to movement can promote the shift back to high prey densities, in what we call an hysteresis break due to movement. In terms of pests control, our results suggest that in order to keep the prey at low densities, predation efforts should be kept as close as possible to $b_{c_2}$.

From a different perspective, if the prey is an endangered specie to be conserved, the regime shift and hysteresis described above may be undesirable since a regime shift can be hard to revert. From the species conservation point of view, it is important to detect regime shifts prior to its occurrence and there has been considerable effort of modelers to find hints that help identifying that the system is close to a regime shift and designing strategies of control to preventing catastrophes [28,34,35,36]. In the present study, we observed that for the upper prey equilibrium $E_2$ (see regions $R_3$ and $R_6$ in Figure 7 and Figure 12, respectively), pattern formation is only observed close to the catastrophic transition value. The typical patterns observed in this case are cold spots what suggests that when the prey species under predation of a generalist predator, starts being scarce in some geographical regions, it may be a signal that some dramatic change is approaching if predation pressure increases. Since the avoidance of regime shifts demands the development and implementation of strategies of control, the indication of a regime shift close to its occurrence may, unfortunately, represent another difficulty of impending it. Our results, thus suggest that monitoring ecological systems is essential for impeding regime shifts.

We hope that our results can, from one hand side, address part of the questions raised by Erbach and colleagues and, on the other side, instigate researchers interested in population dynamics to further investigate generalist predator-prey interactions.