Pattern formation induced by intraspecific interactions in a predator-prey system

TL;DR

This paper investigates how intraspecific interactions like cooperation and competition influence pattern formation in predator-prey systems, revealing new instability mechanisms beyond traditional models with Holling-type functional responses.

Contribution

It introduces a generalized population model with Verhulst-type saturation terms, demonstrating that intraspecific interactions can induce diffusion-driven instability and spatial patterns.

Findings

Intraspecific saturation can cause population instability.

Both cooperative and competitive interactions lead to pattern formation.

Numerical simulations confirm the emergence of spatial patterns.

Abstract

Differential diffusion is a source of instability in population dynamics systems when species diffuse with different rates. Predator-prey systems show this instability only under certain specific conditions, usually requiring Holling-type functionals involved. Here we study the effects of intraspecific cooperation and competition on diffusion-driven instability in a predator-prey system with a different structure. We conduct the analysis on a generalized population dynamics that bounds intraspecific and interspecific interactions with Verhulst-type saturation terms instead of Holling-type functionals. We find that instability occurs due to the intraspecific saturation or intraspecific interactions, both cooperative and competitive. We present numerical simulations and show spatial patterns due to diffusion.

Tables

Table 1. Table 1 : Numerical values used in the simulations shown in Figure 1 , which corresponds to the autocatalytic prey scenario without intraspecific interactions.

Parameters	Numerical values
$r_1$	0.1
$r_2$	0.01
$b_{11}$	0.0
$b_{12}$	-0.00101
$b_{21}$	0.015
$b_{22}$	0.0
$a_1$	0.00001
$a_2$	0.005
$c_1$	0.002
$c_2$	0.005

Table 2. Table 2 : Numerical values used in the simulations shown in Figure  2 and Figure  5 , which correspond to the autocatalytic prey scenario with intraspecific interactions.

Parameters	Numerical values
$r_1$	0.0001
$r_2$	0.6
$b_{11}$	0.0019
$b_{12}$	-0.00075
$b_{21}$	0.00091
$b_{22}$	-0.0019
$a_1$	0.0005
$a_2$	0.000625
$c_1$	0.001251
$c_2$	0.001

Table 3. Table 3 : Numerical values used in the simulations shown in Figure  3 , which corresponds to the autocatalytic prey scenario with intraspecific interactions.

Parameters	Numerical values
$r_1$	0.9
$r_2$	0.00001
$b_{11}$	0.00155
$b_{12}$	-0.001
$b_{21}$	0.00075
$b_{22}$	0.001
$a_1$	0.001
$a_2$	0.001
$c_1$	0.0001
$c_2$	0.0001

Table 4. Table 4 : Numerical values used in the simulations shown in Figure  4 which corresponds to the autocatalytic predators scenario. In Figures  8 - 9 we change b 11 = 0.001915 subscript 𝑏 11 0.001915 b_{11}=0.001915 and in Figures  10 - 11 we used b 11 = 0.001911 subscript 𝑏 11 0.001911 b_{11}=0.001911 .

Parameters	Numerical values
$r_1$	0.02999
$r_2$	-0.090151
$b_{11}$	0.00191
$b_{12}$	-0.0023515
$b_{21}$	0.00105
$b_{22}$	0.0055
$a_1$	0.0021
$a_2$	0.0005
$c_1$	0.001
$c_2$	0.0005

Equations

$$\frac{\partial X_1}{\partial t}$$

$$\frac{\partial X_2}{\partial t}$$

$$f_{11}$$

$$f_{11}{f}_{22}$$

$$f_1\left(u_1,u_2\right)$$

$$f_2\left(u_1,u_2\right)$$

$$f_{11}$$

$$f_{12}$$

$$f_{21}$$

$$f_{22}$$

$$\frac{\max \left(1,q_1\right)}{{\bar{u}}_2}<\mid p_{12}\mid .$$

$$p_{11}{\bar{u}}_1+p_{12}{\bar{u}}_2=\frac{q_1{\bar{u}}_1-1}{1-{\bar{u}}_1},$$

$$\frac{q_1-1}{{\left(1-{\bar{u}}_1\right)}^2}+p_{11}<0.$$

$$p_{21}{\bar{u}}_1+p_{22}{\bar{u}}_2=\frac{q_2{\bar{u}}_2-s}{1-{\bar{u}}_2},$$

$$\frac{q_2-s}{{\left(1-{\bar{u}}_2\right)}^2}+p_{22}>0,$$

$$\frac{q_1-1}{{\left(1-{\bar{u}}_1\right)}^2}+p_{11}>0.$$

$$\frac{q_2-s}{{\left(1-{\bar{u}}_2\right)}^2}+p_{22}<0.$$

$$\delta f_{11}$$

$${\left(\delta f_{11}+f_{22}\right)}^2$$

$${\delta }_c^2{f}_{11}^2+2{\delta }_c\left(2f_{12}{f}_{21}-f_{11}{f}_{22}\right)+f_{22}^2=0.$$

$$p_{21}$$

$$\mid p_{21}{p}_{12}\mid$$

$$\displaystyle\delta_{c}=\biggl{(}\sqrt{\frac{p_{12}}{p_{22}}}$$

$$\frac{\partial X_i}{\partial t}$$

$$z^{\ast }$$

$$\delta$$

$$u_i$$

Full text

Turing pattern \sepreaction-diffusion \seppopulation dynamics \seppredator-prey system \sepcooperation \sepcompetition

††preprint: APS/123-QED

I Introduction

Population ecology treats the increase and fluctuations of populations. Therefore, the purpose of these models is the quantification of the population size of the interact species. In this way, the very first works of Lotka-Volterra equations studied the predator-prey and competing species relations. However, in many of these studies, spatial variation is not considered, but it is necessary element to understand the complete ecological behavior [1].

In particular, Turing instabilities on population dynamics has been studied thoroughly. Many authors have shown that only ecological interactions of opposite sign among species, like predator-prey or parasitism, may produce diffusion driven instability, but pure mutualism or antagonism, with the same sign in interaction between species, may not [2, 1]. Although, a single Lotka-Volterra system can not generate diffusion-driven instability, modified models might. [3] showed that quadratic interactions among populations are needed in order to generate Turing instability in a predator-prey system. They introduced a quadratic positive term for the prey, understood as cooperation and a quadratic negative term for the predators, interpreted as a density dependent death term. Notably, it was also shown that cooperation among predators, introduced as a quadratic expansion of the interaction term, might not produce the same effect. The authors concluded that diffusion-driven instability is caused, in predator-prey system, by self-reinforcement mechanisms acting on the prey, the destabilizers, and self-weaking mechanisms acting on the predators, the stabilizers. [4] also shown that Turing instability might be produced by interference among predators in an innovative way, by constructing a ratio-dependent functional response, using a DeAngelis modified model [5]. [6] and [7] presented another case, using a modified Bazykin model [5], where diffusion-driven instability is also produced by an interference term between predators. In this case, the interference is again produced by a quadratic negative term reflecting predators interference. The authors introduced a prey dependent interaction term between species, instead of a ratio-dependent term. These results somehow contradicts what [8] showed about only ratio-dependent functionals being able of producing diffusion-driven instabilities. Ultimately, [9] showed that using a quadratic term in a Holling-type II functional response also might generate Turing instabilities.

In this paper, we show that another mechanism for Turing instabilities is possible within a predator-prey system. We use a modified version of [10] model to show that using only quadratic interaction terms, adequately bounded by Verhulst-type saturations, may produce diffusion-driven instability. These instabilities appear whether intraspecific direct interactions are allowed or not. When intraspecific direct interactions are not present, the instability arises from the intraspecific saturation acting on the interspecific interaction. When intraspecific direct interactions are allowed, both cooperation and competition terms between predators and preys, promote the instability. All these conditions give rise to different scenarios that we explore in the following section.

II The Model

Diffusion-driven instability takes place in predator-prey systems only under special conditions upon the intraspecific coefficients [1]. For a generic reaction-diffusion system, in dimensionless form, such as:

[TABLE]

it is required, according to [2], that at least the partial derivatives satisfy

[TABLE]

with $f_{ij}=\partial f_{i}/\partial X_{j}$. Here $t$ corresponds to time, the operator $\nabla^{2}$ indicates the Laplacian, the functions $X_{i}$ are the dimensionless populations of the species $i$ and the parameter $\delta$ describes the ratio between their diffusivities ($\delta=d_{2}/d_{1}$). Models with Holling-type II functionals can meet Eqs. (3)–(4) requirements, but Verhulst-type functionals can not meet them [1].

We use a generalized model of population dynamics, based on a modified version of García-Algarra et al. population dynamics model [10, 11], which bounds mutualistic behavior (otherwise unlimited) by saturation Verhulst-like terms. The functionals of a two species system are described with the following equations, in dimensionless form (see A):

[TABLE]

We set $\gamma=1$ to simplify the notation. Let us note that these equations include intraspecific saturation terms, on the environment ($-q_{1}r_{1}$) but also on the interspecific interactions ($1-u_{1}$). The system also allows the existence of $p_{ii}$, which represent direct intraspecific interactions, such as cooperation or competition, which are usually neglected. It is the presence of all these intraspecific terms what allows diffusion-driven instability in a Verhulst-type predator-prey system.

Calculating $f_{ij}$ for the stationary solutions $\bar{u}_{i}$, we obtain:

[TABLE]

Being in dimensionless equations, populations are restricted because of the scaling, within their carrying capacities, to $u_{i}\leq 1$. Without losing generality, we set $u_{1}$ as the prey and $u_{2}$ as the predators from now on. Thus, $p_{12}<0$ and $p_{21}>0$ which mean that $f_{12}<0$ and $f_{21}>0$.

For having diffusion-driven instability, we have two possibilities according to [2]. We might have $f_{11}>0$ and $f_{22}<0$ and we denote this first scenario as the autocatalytic prey. On the other hand, we also might have $f_{11}<0$ and $f_{22}>0$ and we denote this second scenario as the autocatalytic predators. Since the autocatalytic population must be the one which diffuses slower, we have that Eqs. (1)–(2) is coherent with the first scenario, where $d_{2}>d_{1}$. For the second scenario, where $d_{2}<d_{1}$, instead of having $\delta\in\left]1,\infty\right[$ for Turing instability, we have $\delta\in\left]0,1\right[$.

In the absence of terms $p_{ii}$, evaluating Eqs. (5)–(6) for the stationary solutions $\bar{u}_{i}$ force that $\left(1+p_{12}\bar{u}_{2}\right)=\bar{u}_{1}\left(q_{1}+p_{12}\bar{u}_{2}\right)$ and also $\left(s+p_{21}\bar{u}_{1}\right)=\bar{u}_{2}\left(q_{2}+p_{21}\bar{u}_{1}\right)$. Since we already have that $p_{12}<0$, $f_{11}>0$ only occurs if $\mid p_{12}\mid>1/\bar{u}_{2}$ and $\mid p_{12}\mid>q_{1}/\bar{u}_{2}$. This opens the possibility to a new mechanism for diffusion-driven instabilities motivated entirely by the intraspecific saturation on the interspecific interaction of Eqs. (5)–(6). This mechanism corresponds to the autocatalytic prey scenario. No autocatalytic predators scenario is possible, since $u_{2}\left(q_{2}+p_{21}\bar{u}_{1}\right)>0$ for any $q_{2}$ and $p_{21}$.

In the absence of intraspecific saturation on any interactions, either intraspecific or interspecific, Eqs. (7)–(10) give the familiar result of both $f_{11},f_{22}\leq 0$, that does not allow diffusion-driven instability. The case without any intraspecific saturation, not even with the environment, $f_{11},f_{22}=0$, which correspond to the classical Lotka-Volterra system [1].

II.1 Autocatalytic prey without intraspecific interactions

In the absence of intraspecific interactions, i.e. for $p_{ii}=0$, we already saw that $f_{11}>0$ only if $\mid p_{12}\mid>1/\bar{u}_{2}$ and $\mid p_{12}\mid>q_{1}/\bar{u}_{2}$. For this relations, we might derive,

[TABLE]

On the other hand, $f_{22}<0$ always, since $\left(s+p_{21}\bar{u}_{1}\right)>0$ for any $p_{21}$, so no further conditions are needed.

II.2 Autocatalytic prey with intraspecific interactions

In this scenario, and by allowing the presence of $p_{ii}$, we might have $f_{11}>0$ only when $\left(1+p_{12}\bar{u}_{2}+p_{11}\bar{u}_{1}^{2}\right)<0$. From Eqs. (5)–(6) evaluated for $\bar{u}_{i}$, we have,

[TABLE]

so, for $f_{11}<0$ we need that,

[TABLE]

This condition allows two possible behaviors for $p_{11}$, i.e. for the intraspecific interactions of the prey. For $q_{1}>1$, prey must be competitive and $\mid p_{11}\mid>(q_{1}-1)/(1-\bar{u}_{1})^{2}$. For $q_{1}<1$, prey might be competitive, without any restriction, or cooperative, as long as $\mid p_{11}\mid<(q_{1}-1)/(1-\bar{u}_{1})^{2}$. The value of $q_{1}$ comes from their interpretation in Eqs. (1)–(2), through the transformations showed in A. $q_{1}=\left(1/c_{1}\right)/\left(r_{1}/a_{1}\right)$, which might be understood as the ratio between the population limit due exclusively to the resources obtained from the interspecific and intraspecific interactions, $1/c_{1}$, and the population limit due exclusively to the resources from the environment, $r_{1}/a_{1}$.

On the other hand, $f_{22}<0$ requires $\left(s+p_{21}\bar{u}_{1}+p_{22}\bar{u}_{2}^{2}\right)>0$. For the stationary solution, $f_{2}=0$ in Eqs. (5)–(6), we derive,

[TABLE]

so, for $f_{22}>0$ we need that,

[TABLE]

which allows both cooperative and competitive predators, regardless of the sign of $s$. As long as $0<q_{2}-s$, predators might be cooperative, without any restriction, or they might be competitive, as long as $\mid p_{22}\mid<(q_{2}-s)/(1-\bar{u}_{2})^{2}$. But, if $q_{2}-s<0$, predators must be cooperative and $\mid p_{22}\mid>(q_{2}-s)/(1-\bar{u}_{2})^{2}$. We will see later that the intensity of this self-interaction $p_{22}$ will condition the value of the critical diffusion.

II.3 Autocatalytic predators

In this other scenario, the only change is that now $f_{11}<0$ and $f_{22}>0$ are required. For the first condition, it is needed that $\left(1+p_{12}\bar{u}_{2}+p_{11}\bar{u}_{1}^{2}\right)>0$. Using what we derived on the previous section, since Eq. (12) is fulfilled again, instead of Eq. (13), we have,

[TABLE]

Now, this condition allows both cooperative and competitive prey in the opposite direction of what happened in the previous scenario. For $q_{1}>1$, prey might be competitive, as long as $\mid p_{11}\mid<(q_{1}-1)/(1-\bar{u}_{1})^{2}$ and cooperative without any restriction. On the other hand, if $q_{1}<1$, prey must be cooperative and $p_{11}>\mid q_{1}-1\mid(1-\bar{u}_{1})^{2}$. We will also see that the intensity of $p_{11}$ will determine the value of the critical diffusion.

On the other hand, for $f_{22}>0$, we also use Eq. (14) and now, instead of Eq. (15), we have,

[TABLE]

For $q_{2}-s>0$, $p_{22}<0$ and predators must be competitive, but additionally, $\mid p_{22}\mid>(q_{2}-s)/(1-\bar{u}_{2})^{2}$. But, when $q_{2}-s<0$, predator might be competitive, without any restriction, or cooperative, as long as $p_{22}<\mid q_{2}-s\mid/(1-\bar{u}_{2})^{2}$

All the conditions derived in the last two scenarios are only established to see the possible ecological regimes, i.e. the signs allowed on $p_{ii}$, that are valid in order to produce diffusion-driven instability. Since $\bar{u}_{i}$ are functions of $p_{ii}$, no simple relation can be obtained from Eqs. (12)–(17). The same applies to the first scenario, about $p_{12}$ and its relation with $q_{1}$ and $\bar{u}_{2}$. This can be seen in the Results.

II.4 Diffusion-driven instability

Diffusion-driven instabilities require, according to [2], that conditions (3)–(4) change into,

[TABLE]

These conditions explain the reason why, considering the autocatalytic prey scenario with intraspecific interactions, in Eqs. (15) $p_{22}$ was an indicator of the critical diffusion $\delta_{c}$. This applies also for the autocatalytic predators scenario, except that in that case, it is $p_{11}$ the parameter that acts as an indicator of the $\delta_{c}$. From [2] it is required that $\delta_{c}$ follows,

[TABLE]

This means that, at least, $\mid f_{22}/f_{11}\mid<\delta_{c}$ for the autocatalytic prey scenario. In other words, given a $f_{11}$, the greater the cooperation of the predators, the greater critical diffusion will be needed to get a diffusion-driven instability. In the autocatalytic predators scenario we have that $\mid f_{22}/f_{11}\mid>\delta_{c}$, i.e. the inverse dependence is needed between them in order to get the critical diffusion.

We test for diffusion-driven instability using the non-dimensionless system. We introduce small perturbations to the homogeneous stationary solutions of the system, given by $\bar{X}_{i}$. Perturbations are introduced, as functions of fixed wavelength of the form $X_{i}=\bar{X}_{i}+X^{\prime}_{i}e^{\lambda t}e^{ikz}$, into Eqs. (1)–(2) and neglecting the non-linear terms [12]. This gives a set of two equations relating the eigenvalues $\lambda$ with the wavenumber $k$. This constitutes a dispersion relation from which the stability of the system can be verified. We present this relation on the following section, along with the numerical solutions of the nonlinear system.

III Results

III.1 Linear stability analysis

We tested the stability of the homogeneous stationary solution by replacing some test values for $c_{1}$ and $b_{12}$ in the autocatalytic prey scenario without intraspecific interactions, and some test values for $s$ and $p_{ij}$ in both scenarios with intraspecific interactions. Eigenvalues $\lambda$ were obtained as function of wavenumber $k$. $Re\left(\lambda\right)$ change from negative to positive for a certain values of $k$, indicating the cases where a small perturbation with wavelength $2\pi/k$ will not vanish. Instead, those perturbations will grow and will make the system unstable on a linear approximation; the system will stabilize itself by the nonlinear terms.

For the scenario of the autocatalytic prey without intraspecific interactions, we tested the case where both species have a positive dependence on the environment, i.e. $s>0$. In Figure 1, we show the effects of parameters $c_{1}$ and $b_{12}$ in the instability of the system. We see, on the left, that instability is promoted with lower values of $\mid b_{12}\mid$, i.e. with less effects on the prey by the predators and, on the right, with greater values of $c_{1}$, i.e. with a higher intraspecific saturation. However, this effect reaches a point where the system may become intrinsically unstable (see curve C on the right) and no diffusion-driven instability might be generated.

To explore the scenario of autocatalytic prey with intraspecific interactions, we tested the case where both species have a positive dependence on the environment, i.e. $s>0$, and where predators compete and cooperate among themselves. In Figure 2, we show the effects of both competition of predators and cooperation of prey in the instability of the system. We see, on the left, that lower competition among predators promotes a greater instability in the system. But, on the right, we also see that lower cooperation of prey promotes also a greater instability. In Figure 3, we show the effects of both cooperation of predators and cooperation of prey in the instability of the system. On the left, we see the influence of cooperation in prey and how it promotes the system instability, while on the right, we see the influence of cooperation in predators, and how it promotes the stability instead.

For the autocatalytic predators scenario, we tested the case where $s<0$, which means that prey have a positive dependence on the environment, but the predators do not. Also, besides cooperative predators, we use cooperative prey. In Figure 4, we show the effects of both cooperations in the instability of the system. Higher cooperation in both populations promotes Turing pattern formation, but lower cooperation also allows the system to become unstable. Continuing to lower the cooperation further results in a steady state that is unstable even without diffusion.

All solutions we tested were pairs corresponding to saddle-node bifurcations, because they converge and disappear when parameters change [13]. It is interesting to see also that, when parameters change the other way around and solutions diverge, the diffusion-driven instability is lost. Figure 5 shows the phase space of the autocatalytic prey scenario and values correspond to those of Table LABEL:table:autoprey.

The scenarios discussed here reflect some differences with other previous attempts to find diffusion-driven instability in predator-prey systems with cooperative prey. [14] discussed a predator-prey model without saturations and with only the prey depending on the environment. For them, prey were cooperative and predators competitive. Their model required, for diffusion-driven instability to occur, that,

[TABLE]

in the specific scenario of autocatalytic prey with intraspecific interactions. This conditions are not met neither with values in Table LABEL:table:autoprey or Table LABEL:table:autoprey2 .

III.2 Numerical simulations

We solve the nonlinear system by carrying out a numerical simulation of Eqs. (1)–(2). Since only two possible patterns may arise in a one dimensional system, which are identical or inverse [2], we use values of Table LABEL:table:autoprey2 and Table LABEL:table:autopredators to test both scenarios with intraspecific interactions, the autocatalytic prey and the autocatalytic predators. We chose periodic boundary conditions along a one-dimensional space with cell width of size $\Delta z=0.1$ spatial units, where both species $X_{i}$ evolve. We use a simple Euler method with a time step of $\Delta t=0.0001$, which we tested to be accurate. Initial conditions where set with small random perturbations around the homogeneous stationary solutions $\bar{X}_{i}$. Computations where carried out for enough time in order to reach a steady pattern.

For the autocatalytic prey scenario with parameter values of Table LABEL:table:autoprey2, diffusion-driven instability appears with $\delta=45$. The corresponding wavelength of the fastest growth for this diffusion is $29.24$ spatial units. As we use a grid of cells with $\Delta z=0.1$ spatial units, it is expected to have a pattern of three or four peaks in a length of $120$ spatial units. Although an almost uniform pattern of three peaks form for both populations, their amplitudes reduce constantly until reaching a fixed value of $1.35\times 10^{-12}$ for $u_{1}$ and $3.65\times 10^{-13}$ for $u_{2}$. Both dimensionless populations show the same pattern, given that in an autocatalytic prey scenario, both species follow the same dynamics. This is shown in Figure 6. We also show the time evolution of the pattern in Figure 7.

For the autocatalytic predators scenario we conducted two different tests. First, we used parameter values of Table LABEL:table:autopredators and $b_{11}=0.001915$. Diffusion-driven instability appears with $\delta=0.474$. The corresponding wavelength of the fastest growth for this diffusion is $13.49$ spatial units. With $\Delta z=0.1$ spatial units as the cell width, it is expected to have a pattern of nine peaks in a length of $120$ spatial units. Here, the dimensionless populations show an inverse pattern, given that in an autocatalytic predators scenario, species follow the opposite dynamics. This is shown in Figure 8. The spatial pattern is formed with the corresponding wavelength of the Turing instability, but its amplitude continue growing indefinitely. We show the pattern at two different times in Figure 8 and the time evolution in Figure 9.

We also test the case when $b_{11}=0.001911$ and the other parameters where those of Table LABEL:table:autopredators. Diffusion-driven instability appears with $\delta=0.392$ and the corresponding wavelength of the fastest growth for this diffusion is $18.21$ spatial units. For a length of $120$ spatial units, seven peaks would be expected, but we obtained a pattern with nine. This pattern is unstable, and not only its amplitude grows, as in Figure 10 but it oscillates around the stationary solution, as it can be seen in the first steps of Figure 11. Eventually, the amplitudes reaches another basin and the populations go to another stationary solution, the partial extinction of $u_{1}$.

IV Conclusions

Here we studied the generation of patterns from intraspecific interactions, which are usually neglected in most ecological models or which are introduced ad hoc to study specific cases. [15] observed that, among animal species, intraspecific direct interactions act as inhibitory or autocatalytic mechanisms. When individuals behave aggressively among them, this behavior promotes their dispersion across the available territory. On the other hand, when the same individuals cooperate, gregarious behaviors appeared. When these species are involved in an ecological system, both mechanism might couple and diffusion-driven instabilities arise.

In this work, we have shown that intraspecific interactions in a predator-prey system might lead to diffusion-driven instabilities. These intraspecific interactions can be positive (cooperation) or negative (competition), they can act on the predators or on the prey, or even they can be direct (being an active interaction) or indirect (acting as a saturation). This means that they are not as limited as some previous studies pointed out [3, 4, 6, 7, 8, 9]. In the absence of intraspecific direct interactions (terms $b_{ii}X_{i}$), saturation acting on the prey relation with environment resources might cause instability driven by diffusion as long as Eq. (11) holds. This mechanism leads to an autocatalytic prey scenario. No such mechanism exists for predators. When intraspecific direct interactions are present, Turing patterns might arise either with autocatalytic prey or autocatalytic predators, with conditions Eqs. (12)–(17) that allow them to be cooperative or competitive, regardless of the scenario.

We have shown with numerical simulations that instabilities give rise to spatial patterns that might be identical for both species, in the autocatalytic prey scenario, such as Figure 6, or inverse, in the autocatalytic predators scenario, such as Figure 8 or Figure 10. Spatial patterns are only stable in the autocatalytic prey scenario, reaching a fixed amplitude lower than the originally reached. For the autocatalytic predators scenario, initial random perturbations grow continuously, or they show oscillatory patterns of growing amplitude around the stationary solution. Their amplitudes grow until populations reach the basin of an stable stationary solution. Although unstable Turing patterns are known, specially around Hopf bifurcations [16, 17], here we found them around a saddle-node bifurcation. [18] found that Turing instabilities around a saddle-node bifurcation led to stationary Turing patterns, the opposite of what we have found.

Since ecological models deal with living species that inhabit spatial domains, patterns arising from diffusion mechanisms are relevant to a better understanding of the behavior of populations. As [15] pointed out, intraspecific behavior is a well known mechanism to disperse or concentrate individuals. Therefore, pattern formation due to intraspecific interactions should be considered as part of the ecological behaviors that species show among their interactions.

Acknowledgements

This work was supported by Ministry of Economy and Competitiveness of Spain (research project MTM2015-63914-P).

Appendix A Nondimensionalization of the population dynamics system

In this work, we added a diffusion term to a generalized version of the population dynamics model of [10], denoted as,

[TABLE]

We used the following transformation,

[TABLE]

and dropped the ∗ in order to get Eqs. (5)-(6).

Appendix B Numerical values

Here we present the numerical values used in the simulations. We presented the values according to the population equations instead of the dimensionless system, since the latter can be derived from the transformation described in A.