Abstract
In this paper, we consider spatial predator–prey models with diffusion and prey-taxis. We investigate necessary conditions for pattern formation using a variety of non-linear functional responses, linear and non-linear predator death terms, linear and non-linear prey-taxis sensitivities, and logistic growth or growth with an Allee effect for the prey. We identify combinations of the above non-linearities that lead to spatial pattern formation and we give numerical examples. It turns out that prey-taxis stabilizes the system and for large prey-taxis sensitivity we do not observe pattern formation. We also study and find necessary conditions for global stability for a type I functional response, logistic growth for the prey, non-linear predator death terms, and non-linear prey-taxis sensitivity.
1. Introduction
There are basically three mechanisms for spatial pattern formation in systems of two reaction–advection–diffusion equations; the Turing patterns Citation25, chemotaxis patterns Citation12, and patterns created through reaction kinetics, e.g. the Brusselator Citation15. Turing patterns typically arise for a fast inhibitor and a slow activator. Chemotaxis patterns are based upon aggregation towards a chemical signal.
For a predator–prey system without prey-taxis, Okubo and Levin Citation26 note that an Allee effect in the functional response and a density-dependent death rate of the predator are necessary to generate spatial patterns. The inclusion of species migration (constant flow) as an additional transport process may also increase the possibility of pattern formation Citation13. The directional movement of zooplankton plays a role in generating patterns in a plankton community model Citation22. For the same system, diffusion also generates pattern formation, and the combined effects of diffusion and velocity result in spatial pattern and an instability like travelling waves, i.e. travelling patchy distributions Citation23. Indeed, the magnitude of the relative flow velocity determines the flow-induced instability Citation28. Various travelling wave solutions have been studied with similar systems Citation4 Citation18 Citation33. In particular, Citation4 Citation33 also considered escaping behaviours of prey from predation.
For chemotaxis models, spatial patterns have been studied analytically and numerically (see, for example, Citation7 Citation35 Citation36). In contrast to the rich development of chemotaxis models, the pattern formation of prey-taxis models is still open to wide investigations. Lewis Citation20 studied pattern formation in plant and herbivore dynamics and herbivory-taxis was seen to reduce the likelihood of pattern formation. Arditi et al. Citation3 and Chakraborty et al. Citation6 considered a different aspect of predator response to the prey distribution that the velocity of the predators is dependent on prey density.
The goal of this paper is to investigate the contribution of predator and prey movements to spatial pattern formation in predator–prey systems. In particular, we consider foraging behaviour of predators that move towards high prey density. For that, we extend the predator–prey diffusion–reaction model in Citation27 by incorporating the concept of prey-taxis Citation14.
1.1. The model
A prey-taxis model was derived by Kareiva and Odell in Citation14, and they studied predator aggregation in high prey density areas. Later the model was applied to estimate the mean travel time of a predator to reach a prey resource Citation9. Here we extend the Kareiva and Odell model to studying pattern formation.
The prey-taxis model discussed below contains both diffusion terms that might lead to Turing-type instabilities and a prey-taxis term that might lead to aggregation of predators on local concentrations of prey. In this paper, we will investigate the relative importance of these effects for spatial pattern formation. Prey-taxis allows predators to search more actively for prey, and can generate different spatial patterns from those formed in models without prey-taxis. Generally speaking, we find that prey-taxis tends to stabilize the predator–prey interactions.
The characteristic feature of prey-taxis equations is that taxis is incorporated as an advection term Citation14 Citation17. In this paper, we consider the following prey-taxis model
To investigate the pattern formation properties of EquationEquations (1) and Equation(2)
, we first consider Equation(1)
and Equation(2)
without taxis, i.e. χ=0. Secondly, we study the full model Equation(1)
and Equation(2)
with χ≠0. In Section 2.1, we study pattern formation for EquationEquations (1)
and Equation(2)
. It turns out that pattern formation crucially depends on the functional forms of functional response h(v, n), on the death rate δ (n), and on the prey growth kinetics f(v). We investigate typical cases, that are discussed in the literature, see, e.g. Citation34. For h(v, n) we consider type I (linear) functional response h(v, n)=v, type II (hyperbolic) functional response
, linear-ratio functional response
, and hyperbolic-ratio functional response
. The death rate δ(n) is either constant
or density-dependent
. For the prey kinetics, we assume either logistic growth f(v)=1−v or an Allee effect
. The above parameters α, ν0, μ, d, ν, K, and a are all positive constants. We summarize the choices of these functions, the corresponding pattern formation results, and the corresponding section in .
Table 1. The possibility of spatial pattern formation is considered in the spatial predator–prey system Equation(1)
and Equation(2)
with various functional responses, h, prey population dynamics, f, and predator death rates, δ.
In Section 3 we consider global stability of the system Equation(1) and Equation(2)
with a type I functional response, density-dependent predator death rate, logistic prey growth rate, and a prey-taxis term. We construct a Lyapunov functional and find that for some condition the coexistence steady state is globally stable. We finish the paper with a discussion and suggestions for further studies (Section 4).
Note that in this paper we implement efficient and accurate numerical methods for each term via a fractional step method Citation19 Citation36 by using MATLAB. For diffusion and reactions terms, we use the Crank–Nicolson scheme and a second-order Runge–Kutta scheme, respectively Citation1 Citation32. For the advection term, we use a high-resolution central scheme Citation16.
2. Pattern formation in prey-taxis systems
In this section we focus on constant prey-taxis and study EquationEquations (1)
and Equation(2)
on an interval [0, L] with homogeneous Neumann boundary conditions given by
We first consider EquationEquations (1) and Equation(2)
for general f(v), h(v, n), and δ (n) and study the specific functional forms later. Since we are interested in understanding biological phenomena, the prey growth function f(v) can be negative with an Allee effect, but the functional response h(v, n) is assumed non-negative. We assume that a non-trivial coexistence steady state
exists.
In order to investigate pattern formation, we follow the standard Turing stability analysis (see Citation25 Citation17 for details).
We first assume that is linearly stable for the purely kinetic equations.
Assumption
Now, we consider the full reaction–taxis–diffusion system Equation(1) and Equation(2)
and obtain the following characteristic equation for an eigenvalue λ of the linearization at
:
Lemma 2.1
Assume that A, B, C, and D are defined in Equation
Equation(5)
and ε is positive. In addition, if AD−BC>0 is assumed as in Equation
Equation(4)
, then there exists χ*≥0 such that M
2(k
2)>0 for all k and all
. In this case, the homogeneous solution
is linearly stable.
Proof
Since AD−BC>0 and ε>0, positive guarantees that M
2(k
2)>0. From the expression
, we can isolate χ and set this χ as χ0. Then we have
. We now define
. We have χ*≥0 and for all
, we have M
2(k
2)>0 independent of the value k. ▪
Therefore, prey-taxis tends to reduce the occurrence of dispersal-induced instability. It is indeed the predator diffusion that is crucial to dispersal-induced instability. When prey act anti-predator defensive behaviours such as kicking and attacking and show chemical defences Citation21, predators may retreat from high prey area, in which case χ in EquationEquation (8) can be negative. As a result, a predator–prey system may generate pattern formation. But here we do not consider this case in detail. In the absence of predators, the prey diffusion would reduce local prey maxima and equilibrate the prey distribution. If predators are present and if they are attracted to local prey maxima through prey-taxis, then the reduction of local prey maxima is enhanced. Hence, if the taxis component is strong enough, we might not expect pattern formation. This is indeed the case, as given in Lemma 2.1. For a specific example, we refer to Section 2.2.
In the following subsections, we consider specific choices for the functional responses, h, the death rate of the predator, δ, and prey growth rate, f. The ability of the prey taxis model Equation(1) and Equation(2)
to exhibit a spatial pattern crucially depends on the parameter functions h(v, n), δ (n), and f(v). Thus, in this section, we study various typical cases separately. An overview of the cases and the corresponding results is given in .
2.1. Type I functional response, density-dependent predator death rate, and Allee effect with diffusion only
In this subsection, we show that pattern formation is possible when there is a type I functional response and an Allee effect along with a density-dependent predator death rate (, row 1). We consider an Allee effect on the prey population dynamics with 0<a<1 and K=4/(1−a)2, a type I functional response, h(v, n)=v, a density-dependent predator death rate,
, ν≥0, and no taxis, i.e. χ=0. Here the parameter a is a threshold, below which the prey population declines. Okubo and Levin Citation26 argued that a predator–prey model with dispersal may generate diffusion-driven instability if the mortality of the predator depends on the population density and the per-capita growth rate of the prey is determined by an Allee effect. Note that the trivial steady state (v, n)=(0, 0) is locally stable because for (v, n)=(0, 0) the characteristic polynomial for purely kinetic equations has two negative eigenvalues,
and λ=−K a. We assume biologically relevant parameters in the region 0<δ<1 and 0<a<1. For the prey-only steady state (v, n)=(1, 0), the characteristic polynomial has one positive eigenvalue
and one negative eigenvalue
.
For the homogeneous coexistence steady state , we find
Here we set . Hence, when
, A is negative and when
, A is positive. Recall that the coexistence steady state
comes from the intersection of the two nullclines:
and
(see ).
Figure 1. The v-nullcline N=K(1−V)(V−a) is shown as a solid curve. For two values of ν, we show the corresponding n-nullcline, N=(V−δ)/ν as a dashed line and a solid line. The equilibrium (v s, n s) is the intersection of the nullclines. We have chosen two values of ν so that v s<v¯ for ν1 and v s>v¯ for ν2 with ν1<ν2.
![Figure 1. The v-nullcline N=K(1−V)(V−a) is shown as a solid curve. For two values of ν, we show the corresponding n-nullcline, N=(V−δ)/ν as a dashed line and a solid line. The equilibrium (v s, n s) is the intersection of the nullclines. We have chosen two values of ν so that v s<v¯ for ν1 and v s>v¯ for ν2 with ν1<ν2.](/cms/asset/f7f0fce2-daa7-4976-ad7f-ad5648d1bf17/tjbd_a_371781_o_f0001g.gif)
First, we consider . Since
, we find that at v=v
s the v-nullcline is above the n-nullcline. This means that
, which translates into the condition
For this case we prove stability.
Lemma 2.2
Assume that h(v, n)=v,
and χ=0. If
then no pattern formation occurs about the coexistence steady state
for the system
Equation(1)
and
Equation(2)
.
Proof
The condition implies that
. Hence A<0. In addition, we find B<0, C>0, and D<0 and A<0, B<0, and D<0 imply
for all real k. Hence, we cannot expect diffusion-taxis-driven instability about the coexistence steady state. ▪
In it is noted that v
s should be between a and 1, i.e. a<v
s<1, otherwise n
s is negative. In Lemma 2.2, we considered that and found no pattern. Thus we now consider
. First we investigate how many v
s may exist between a and v¯, and then we find conditions for the existence of v
s between a and v¯.
The v values for the coexistence steady state are obtained from
When a<δ, shows that EquationEquation (12) has two real roots. Indeed, for the root less than a, n
s would be negative, which is not biologically relevant. Hence, when a<δ, EquationEquation (12)
has one biologically relevant root. In addition,
leads to
. Therefore, for
When a>δ, we may expect two positive roots from EquationEquation (12). However, a simple computation of EquationEquation (12)
shows that we cannot have two positive roots. Under assumption Equation(13)
, the biologically relevant solution of EquationEquation (12)
is given by
The discriminant in EquationEquation (14) is zero for
We find A>0 from the condition Equation(13). Additionally, from EquationEquation (9)
we find B<0, C>0, and D<0. The stability condition A+D<0 leads to a condition
In we plot the left- and right-hand sides of EquationEquation (16) as a function of v
s.
Figure 2. Plot of the left- and right-hand sides of EquationEquation (16) as function of v
s. The region of A+D<0 is where the dashed line lies above the curve.
![Figure 2. Plot of the left- and right-hand sides of EquationEquation (16) as function of v s. The region of A+D<0 is where the dashed line lies above the curve.](/cms/asset/ab1190a9-be37-4735-b93d-b4e9ea5fcacd/tjbd_a_371781_o_f0002g.gif)
As γ varies from zero to infinity, the intersection of and
changes from
to
. Given a value for v
s, we can always choose γ small enough such that condition Equation(16)
is not true. Thus, γ should be greater than a minimum value, γ0. Here
, where v
s is computed in EquationEquation (14)
. Therefore, for
, we have A+D<0.
Thus a biologically relevant v s is in the interval
We found that EquationEquation (17) holds under assumption Equation(13)
.
Theorem 2.3
Assume that h(v, n)=v,
and χ=0. If a satisfies condition
Equation(13)
, then (i) the coexistence steady state
exists (ii) AD−BC>0, (iii) if in addition, there exists
such that for each
there exists a non-empty interval [k
1, k
2] of unstable modes, so we may expect diffusion-driven instability about the coexistence steady state, and (iv) if
then
is linearly stable.
Proof
-
(i) It was shown that condition Equation(13)
implies the existence of a unique positive root v s.
-
(ii) When a positive v s exists, v s satisfies condition Equation(17)
. Now we consider the condition for AD−BC>0.
. Indeed, this is true by condition Equation(17)
. Therefore, AD−BC is always positive under the assumption of the existence of a coexistence steady state.
-
(iii) M 1(k 2) and M 2(k 2) are given by EquationEquations (10)
and Equation(11)
, respectively, with
, B=−v s<0,
, and
. Hence, A+D<0 guarantees
. If
, then D<0 gives
, hence M 2(k 2) is always positive. On the other hand, setting
, for
, we have
and M 2(k 2) can be negative for some k. Setting T=k 2, the quadratic equation M 2(T)=0 may have two roots, T 1, 2 (see ). By solving this quadratic equation, it can be shown that for
there exist real k 1 and k 2. For unstable modes k∈[k 1, k 2 with
and
, we have Re(λ)>0. Hence we may expect diffusion-driven instability about the coexistence steady state (see also Citation30).
-
(iv) If
, then M 2(k 2) is always positive for all k. Hence we cannot expect diffusion-driven instability about the coexistence steady state.
Segel and Jackson Citation30 also considered diffusion-driven instability in a predator–prey interaction. They used and f(v)=1+Kv, and found the wavelength of the instability (see also Citation26 for general discussion on diffusion-driven instability in a predator–prey interaction).
In particular, for , we apply a perturbation method to approximate two values T
1, 2, that is,
and
. Therefore, for
, M
2(T) is negative.
For example, we consider an interval [0, L] with homogeneous Neumann boundary condition given by EquationEquation (3). If k(n)=nπ/L in [k
1, k
2] with positive integer n, then pattern formation occurs. Thus, we can calculate a minimum domain size for pattern formation. Since
, we substitute k(n)=nπ/L and rearrange the inequality with respect to L. Then we have
In , we show phase portraits of the predator–prey system Equation(1) and Equation(2)
without dispersal terms. As γ increases, the coexistence steady state bifurcates from an unstable spiral to a stable spiral. From simulations with various γ, it is noted that an unstable limit cycle occurs for a certain range of γ. When γ is smaller than the lower bound of this range, the coexistence steady state is an unstable spiral. When γ is bigger than the upper bound of the range, the coexistence steady state is a stable spiral with non-empty basin of attraction. shows that the stable coexistence steady state without dispersal terms becomes unstable if diffusion terms are introduced. As a result, patterns are generated. We demonstrate a snapshot of the asymptotic prey and predator distributions in . It is noted that high prey density area seems to attract more predators. Moreover, the patch size of prey is seen to be an important factor to attract more predators.
Figure 4. Coexistence steady state is shown to be locally asymptotically stable for system Equation(1) and Equation(2)
without dispersal terms and with f(v)=16(1−v)(v−0.5), h(v, n)=v, and δ (n)=0.6+0.1 n. Time step is d t=0.01 and γ=13. Here the coexistence steady state is (v
s, n
s)=(0.695, 0.952). Available in colour online.
![Figure 4. Coexistence steady state is shown to be locally asymptotically stable for system Equation(1) and Equation(2) without dispersal terms and with f(v)=16(1−v)(v−0.5), h(v, n)=v, and δ (n)=0.6+0.1 n. Time step is d t=0.01 and γ=13. Here the coexistence steady state is (v s, n s)=(0.695, 0.952). Available in colour online.](/cms/asset/52b6821f-d20a-453a-beeb-319ab5bc55c5/tjbd_a_371781_o_f0004g.jpg)
Figure 5. Coexistence steady state is shown to be locally unstable for system Equation(1) and Equation(2)
with χ(v)=0.0, f(v)=16(v−0.5)(1−v), h(v, n)=v, and δ (n)=0.6+0.1 n. Spatial grid size is d x=0.25, time step d t=0.01, and γ=14 with 60 time units. The diffusion coefficient ε is 0.01. Here the coexistence steady state is (v
s, n
s)=(0.695, 0.952). Available in colour online.
![Figure 5. Coexistence steady state is shown to be locally unstable for system Equation(1) and Equation(2) with χ(v)=0.0, f(v)=16(v−0.5)(1−v), h(v, n)=v, and δ (n)=0.6+0.1 n. Spatial grid size is d x=0.25, time step d t=0.01, and γ=14 with 60 time units. The diffusion coefficient ε is 0.01. Here the coexistence steady state is (v s, n s)=(0.695, 0.952). Available in colour online.](/cms/asset/3a14f709-cbd3-4ccb-aebb-d2806a115791/tjbd_a_371781_o_f0005g.jpg)
2.2. Type I functional response, density-dependent predator death rate, and Allee effect with diffusion and prey-taxis
We now include prey-taxis in the calculations of the previous subsection (, row 1). We consider the reaction–diffusion–taxis system Equation(1) and Equation(2)
for
. We consider Allee-type growth for the prey,
with 0<a<1 and K=4/(1−a)2, a type I functional response, h(v, n)=v, and a density-dependent predator death rate,
. We have shown in the previous subsection that for χ=0, pattern formation may occur. In this subsection, we consider how the conditions of pattern formation change if χ is introduced.
Lemma 2.4
Assume a, δ, and ε satisfy instability conditions of Theorem 2.3. Then from Lemma 2.1 we compute
Proof
Here M
1(k
2) from Equation(7) is the same as in the case of diffusion-only Equation(10)
so that it is negative for all k. But M
2(k
2) from Equation(8)
is different by the term Bχ n
s. Setting M
2(k
2)=0 and T=k
2, we obtain after rearrangements
Figure 7. Plot of the left- and right-hand sides of EquationEquation (18) as a function of T with T=k
2. The solid curve is from the left-hand side of EquationEquation (18)
and the dashed lines are from the right-hand side of EquationEquation (18)
. As χ increases, the number of intersection changes from zero to two. Note that B is negative.
![Figure 7. Plot of the left- and right-hand sides of EquationEquation (18) as a function of T with T=k 2. The solid curve is from the left-hand side of EquationEquation (18) and the dashed lines are from the right-hand side of EquationEquation (18). As χ increases, the number of intersection changes from zero to two. Note that B is negative.](/cms/asset/72d249c0-c34d-47bf-b05a-fce8ad0becbe/tjbd_a_371781_o_f0007g.gif)
shows that the coexistence steady state for the spatially homogeneous predator–prey system Equation(1) and Equation(2)
without dispersal terms is stable. In , introducing the diffusion term generates patterns. The numerical simulations confirmed that when we introduce a large prey-taxis term, patterns disappear (not shown here).
2.3. Linear ratio-dependent functional response, constant predator death rate, and logistic growth
In this subsection, we show that pattern formation is impossible when there is a linear ratio functional response and logistic growth along with a constant predator death rate (, row 2). We consider the linear ratio-dependent functional response, with logistic growth for the prey, f(v)=1−v, and a constant predator death rate,
, and ν0 is a constant parameter. Thus the coexistence steady state is now
, which is biologically relevant for
. In this case, we obtain
We observe that A<0, D<0, B=0, AD−BC>0, and for all k. Hence, the homogeneous steady state is linearly stable.
Lemma 2.5
Assume f(v)=1−v,
and
then no pattern formation occurs about the coexistence steady state,
for system
Equation(1)
and
Equation(2)
.
2.4. Hyperbolic ratio-dependent functional response, constant predator death rate, and logistic growth
We now modify the analysis of the previous subsection to include a hyperbolic ratio rather than linear ratio functional response. This allows for the possibility of pattern formation, providing that taxis is sufficiently small (, row 3). We consider hyperbolic ratio-dependent functional response, with logistic growth for the prey, f(v)=1−v, and a constant predator death rate,
, and μ≥0 and d≥0 are constants. Thus the coexistence steady state in this case is
We consider conditions that A+D<0 and AD−BC>0. For A<0, it is seen that A+D<0 and AD−BC>0. For A>0, , with
Lemma 2.6
Assume f(v)=1−v,
and χ=0.
-
(i) If
no pattern formation occurs about the coexistence steady state,
for system Equation(1)
and Equation(2)
.
-
(ii) Assume
. There exists
such that for each
there exists a non-empty interval [k 1, k 2 of unstable modes, so we may expect diffusion-driven instability about the coexistence steady state,
-
(iii) in case (ii) if
then
is linearly stable.
Proof
-
(i) First,
implies A<0. In addition, B<0, C>0, and D<0 result in positive M 2(k 2) from Equation(11)
. Hence, we cannot expect diffusion-taxis-driven instability about the coexistence steady state.
-
(ii) Second, we consider
, which gives A>0. It is also seen that AD−BC>0 and for
, A+D<0, which implies that M 1(k 2) from Equation(10)
is negative. However, when ε is less than
, then A+ε D is positive. Thus, M 2(k 2) can be negative. With the same steps in Theorem 2.3, we can find k 1 and k 2 with
, where
, there exist real k 1 and k 2. Furthermore, for k 1<k<k 2, we have Re(λ)>0, and we may expect diffusion-driven instability about the coexistence steady state.
-
(iii) If
, then M 2(k 2) is positive for all k. Hence, we cannot expect diffusion-driven instability about the coexistence steady state.
Alonso et al. Citation2 also considered a hyperbolic ratio-dependent functional response for pattern formation by using numerical exploration of the parameter space.
Now we can follow the argument of the case including an Allee effect. Thus, the reaction–diffusion system may show diffusion-driven instability depending on parameters μ, d, δ, γ, and ε. Furthermore, the prey-taxis term tends to inhibit the occurrence of dispersal-driven instability (see Lemma 2.1 and Section 2.2 for the full argument).
In , we show phase portraits of the predator–prey system Equation(1) and Equation(2)
with hyperbolic-ratio functional response and without dispersal terms. As γ increases, the coexistence steady state bifurcates from an unstable spiral to a stable spiral. demonstrates that this homogeneous coexistence steady state becomes unstable if diffusion terms are introduced. As a result, patterns are generated. It is shown that when we introduce a large prey-taxis term, patterns eventually disappear (see Citation17 for figure).
Figure 8. Coexistence steady state is shown to be locally asymptotically stable for system Equation(1) and Equation(2)
without dispersal terms and with h(v, n)=0.8 v/(0.05 n+v), f(v)=1−v, and δ (n)=0.76. Time step is d t=0.005, and γ=15. Here the coexistence steady state is (v
s, n
s)=(0.2, 0.211). Available in colour online.
![Figure 8. Coexistence steady state is shown to be locally asymptotically stable for system Equation(1) and Equation(2) without dispersal terms and with h(v, n)=0.8 v/(0.05 n+v), f(v)=1−v, and δ (n)=0.76. Time step is d t=0.005, and γ=15. Here the coexistence steady state is (v s, n s)=(0.2, 0.211). Available in colour online.](/cms/asset/749089d6-8a46-4d2b-8ccb-aa207b4152dc/tjbd_a_371781_o_f0008g.jpg)
Figure 9. Coexistence steady state is shown to be locally unstable for system Equation(1) and Equation(2)
with h(v, n)=0.8 v/(0.05 n+v), f(v)=1−v, and δ (n)=0.76 and with χ(v)=0.0. The diffusion coefficient ε is 0.01. Spatial grid size is d x=0.25, time step dt=0.01, and γ=15. Here the coexistence steady state is (v
s, n
s)=(0.2, 0.211). Available in colour online.
![Figure 9. Coexistence steady state is shown to be locally unstable for system Equation(1) and Equation(2) with h(v, n)=0.8 v/(0.05 n+v), f(v)=1−v, and δ (n)=0.76 and with χ(v)=0.0. The diffusion coefficient ε is 0.01. Spatial grid size is d x=0.25, time step dt=0.01, and γ=15. Here the coexistence steady state is (v s, n s)=(0.2, 0.211). Available in colour online.](/cms/asset/4f49c1fe-4fae-4b07-8a97-3272de2b8f9c/tjbd_a_371781_o_f0009g.jpg)
2.5. Type II functional response, density-dependent predator death rate, and logistic growth
Next we consider a hyperbolic functional response and a density-dependent predator death rate from the setting of the previous subsection and show that pattern formation is possible, provided that taxis is sufficiently small (, row 4). We consider type II functional response, as in Citation27, a density-dependent predator death rate,
, and logistic growth for the prey, f(v)=1−v with α>0, 0<δ<1, and ν>0. The coexistence steady state can be obtained from the root of the following system
By applying the intermediate-value theorem, it is shown that there is at least one point v=v
s in the open interval (0, 1) such that g(v
s)=0. Moreover, since for v
s∈(0, 1), n
s corresponding to v
s is positive as well.
Lemma 2.7
Assume f(v)=1−v,
and
then there exists at least one coexistence steady state,
for the system
Equation(1)
and
Equation(2)
.
For an homogeneous coexistence steady state , we find
It is noted that B<0, C>0, and D<0. Thus for A<0, we cannot expect spatial pattern because M 1(k 2)<0 and M 2(k 2)>0 for all k. The condition for A<0 is rewritten in terms of parameters α, δ, and ν as follows:
Therefore, we can summarize the result.
Lemma 2.8
Assume f(v)=1−v,
and
. If
then no pattern formation occurs about the coexistence steady state,
for the system
Equation(1)
and
Equation(2)
.
This result was also confirmed numerically for selected parameter values (not shown here).
Now we consider the case of A>0, that is,
The condition for
Under condition Equation(19), we set the positive root of A+D expressed above with v*, and then v* is between 0 and
. Thus, for v
s in
, we have A+D<0. Note:
. Recall that v
s is independent of γ, so by controlling γ we can make v* smaller than v
s.
Now we consider AD−BC>0. Because A is positive, we cannot guarantee AD−BC>0. After rearrangement, we find
Theorem 2.9
Assume that
f(v)=1−v, and χ=0. {\rm (i)} Then at least one coexistence steady state
exists; (ii) there exists a
such that for all
we have A+D<0 and AD−BC>0; (iii) there exists
such that for each
there exists a non-empty interval [k
1, k
2] of unstable modes, so we may expect diffusion-driven instability about the coexistence steady state; (iv)} if
then
is linearly stable.
Proof
Property (i) was shown in Lemma 2.7. (ii) It was also shown that for , we have A+D<0 and for
, AD−BC>0. Hence, we define
, so for
, we have A+D<0 and AD−BC>0. (iii) Now M
2(k
2) is
If , then D<0 gives
, hence M
2(k
2) is always positive, which results in no diffusion-driven instability for
. Therefore, ε should be strictly less than 1. Indeed, setting
, then for
, we have
and M
2(k
2) can be negative for some k.
By setting T=k
2, we have a quadratic form of M
2(T), that is, . Solving this quadratic form gives that there are two positive ε, say ε1 and ε2 with
, such that for
, M
2(T) has two positive roots, say T
1 and T
2 with T
1<T
2. Therefore, M
2(T) is negative for k∈(k
1, k
2) with
and
. For unstable modes k∈[k
1, k
2], we have Re(λ)>0. Hence we may expect diffusion-driven instability about a coexistence steady state
.
(iv) If , then M
2(k
2) is always positive for all k. Hence, we cannot expect diffusion-driven instability about the coexistence steady state. ▪
For χ≠0, we obtain a similar result as in Section 2.2.
Lemma 2.10
Assume that instability conditions of Theorem 2.9 are satisfied. Then we define
Proof
See Lemma 2.1 ▪
demonstrates that the stable coexistence steady state without dispersal terms becomes unstable when diffusion terms are introduced (Theorem 2.9). As a result, patterns are generated. It is shown that when we introduce a large prey-taxis term, patterns are disappearing (Lemma 2.10 and see Citation17 for figure).
Figure 10. Coexistence steady state is shown to be asymptotically unstable for system Equation(1) and Equation(2)
with χ=0, f(v)=1−v, h(v, n)=v(α+1)/(α+v), and δ(n)=δ+ν n with α=0.2, δ=0.6, γ=1.2, and ν=0.4. The diffusion coefficient ε is 0.01. Spatial grid size d x=0.25 and time step d t=0.01 with 60 time units. Here the coexistence steady state is (v
s, n
s)=(0.296, 0.291). Available in colour online.
![Figure 10. Coexistence steady state is shown to be asymptotically unstable for system Equation(1) and Equation(2) with χ=0, f(v)=1−v, h(v, n)=v(α+1)/(α+v), and δ(n)=δ+ν n with α=0.2, δ=0.6, γ=1.2, and ν=0.4. The diffusion coefficient ε is 0.01. Spatial grid size d x=0.25 and time step d t=0.01 with 60 time units. Here the coexistence steady state is (v s, n s)=(0.296, 0.291). Available in colour online.](/cms/asset/47827c75-7e21-495b-a7d4-42f1a43cb996/tjbd_a_371781_o_f0010g.jpg)
2.6. Type II functional response, constant predator death rate, and logistic growth
We now modify the analysis of the previous subsection to include a constant rather than density-dependent predator death rate. This does not allow for the possibility of pattern formation (, row 5). We consider type II functional response, as in Citation27, a constant predator death rate,
, and logistic growth for the prey, f(v)=1−v. Thus the coexistence steady state is
It is noted that type II functional response does not play any role in pattern formation versus type I. In a numerical solution (not shown) with χ=6.5 and a randomly chosen initial distribution, we observe that the solution converges to the coexistence equilibrium .
2.7. Type I functional response, density-dependent predator death rate, and logistic growth
We now modify the analysis of Section 2.5 to include a type I rather than type II functional response. This also does not allow for the possibility of pattern formation (, row 6). We include competition in the predator death rate, so the predator death rate is . In addition, we consider type I functional response, h(v, n)=v, and logistic growth for the prey, f(v)=1−v. Thus the coexistence steady state is
, which is biologically relevant for
. In this case, we obtain
We find A<0 and D<0 for biologically relevant δ, which result in A+D<0. Moreover, B<0 and C>0 give rise to AD−BC>0. In addition, A<0, D<0, B<0, and AD−BC>0 give for all k. Hence, we note that M
2(k
2)>0 for all k, hence the homogeneous steady state is linearly stable.
This result was also confirmed numerically for selected parameter values (not shown here).
3. Global stability
In the previous section, we showed that without both the Allee effect and the density-dependent predator death rate, diffusion, and prey-taxis do not change the local stability of the coexistence steady state. We choose one of the cases without pattern formation to study the global stability of . We consider logistic growth for the prey, f(v)=1−v, a type I functional response, h(v, n)=v, a density-dependent predator death rate,
, and a non-constant prey sensitivity, χ (v)=b/v, for the spatially homogeneous case of system Equation(1)
and Equation(2)
on an interval Ω=[0, L] with homogeneous Neumann boundary conditions Equation(3)
. The following Lyapunov function,
Theorem 3.1
Let f(v)=1−v, h(v, n)=v,
and
with boundary condition
Equation(3)
. We assume that
then the functional V(v, n) defined in Equation
Equation(20)
is a strong Lyapunov function for system
Equation(1)
and
Equation(2)
. The sets
are positively invariant and
is globally asymptotically stable.
Proof
Setting for L large enough, then we claim that the sets N
L
are positive invariant. When
,
becomes zero due to [Vtilde]=0. For v>v
s, the first term of EquationEquation (20)
is positive. Similarly, the second term of Equation(20)
is positive. Therefore, [Vtilde] is positive in N
L
. Hence, the functional, V(v, n), is bounded below by zero. Moreover, the definition of the Lyapunov functional,
4. Conclusion
In this paper, we considered pattern formation for a predator–prey taxis model of reaction–diffusion–advection type given by Equation(1) and Equation(2)
. We considered various reaction terms: for the predator term they include type I and type II functional responses as well as ratio-dependent functional responses. We considered constant and density-dependent death rate of the predator and logistic growth or an Allee-type growth for the prey.
In summary, the following functional forms support spatial pattern formation:
-
(i) a density-dependent death rate, e.g.
, and an Allee effect, e.g.
, and a type I functional response, e.g. h(v, n)=v:
-
patterns form with no prey-taxis (Section 2.1),
-
patterns persist with small prey-taxis but disappears for large prey-taxis (Section 2.2),
-
patterns disappear when Allee dynamics are replaced by logistic dynamics (Section 2.7),
-
-
(ii) a hyperbolic ratio-dependent functional response, e.g.
and logistic growth, e.g. f(v)=1−v:
-
patterns form with no prey-taxis (Section 2.4),
-
patterns persist with small prey-taxis but disappears for large prey-taxis (Section 2.4),
-
patterns disappear if a hyperbolic functional response is replaced by a linear functional response (Section 2.3);
-
-
(iii) a density-dependent death rate, e.g.
, and a type II functional response, e.g.
.
-
patterns form with no prey-taxis (Section 2.5),
-
patterns persist with small prey-taxis but disappears for large prey-taxis (Section 2.5),
-
patterns disappear if a density-dependent predator death rate is replaced by a constant predator death rate (Section 2.6).
-
The significance of this research is as follows; contrary to a diffusion process that may give rise to pattern formation, prey-taxis tends to stabilize predator–prey interactions (Theorem 2.3). In the long run, prey-taxis tends to transform heterogeneous environments into homogeneous environments, which gives an opposite result to the chemotaxis case. Under strong chemotactic sensitivity, amoebae tend to aggregate Citation29. Hence the role of taxis may be strongly related to the local population dynamics of the species.
In this paper, prey-taxis is shown to tend to reduce the likelihood of pattern formation in spatial predator–prey systems, but other kinds of taxis may have the opposite effect on pattern formation. For example, we may investigate prey defences. Prey tend to adjust their relative position to the predator to reduce predation risk Citation10
Citation24
Citation37
Citation38. We may apply the concept of prey-taxis to prey escape response to predator density. It may refer to predator-taxis. For instance, crayfish (prey) exhibit different activities depending on the presence of a predator (bass). An increased predation risk restricts crayfish foraging and increases anti-predator behaviour such as shelter seeking Citation8
Citation11. Another interesting taxis is that predators may attract their prey to come nearby Citation31. In this case, prey move towards predators. As a conjecture, from EquationEquation (8), we may predict that positive predator-taxis (away from predators) tends to generate pattern formation but negative predator-taxis (towards predators) tends to inhibit pattern formation. However, the detailed argument is left for future work.
Acknowledgements
We thank G. de Vries, C. Cosner, and C.J. Bampfylde for helpful discussion and comments on the manuscript. We appreciate two anonymous reviewers for helpful comments and suggestions. We would also like to thank all the members of the Lewis lab for encouragement and feedback, which have improved this research. JML was supported by a University of Alberta Studentship and by an NSERC collaborative Research Opportunity Grant. TH was supported by NSERC. ML was supported by NSERC and a Canada Research Chair.
References
- Allen , M. B. and Eli , M. B. 1998 . Numerical Analysis for Applied Science , New York : Wiley-Interscience .
- Alonso , D. , Bartemeus , F. and Catalan , J. 2002 . Mutual interference between predators can give rise to Turing spatial patterns . Ecology , 83 : 28 – 34 .
- Arditi , R. , Tyutyunov , Y. , Morgulis , A. and Govorukhin , V. 2001 . Directed movement of predators and the emergence of density-dependence in predator-prey models . Theoret. Population Biol. , 59 : 207 – 221 .
- Biktashev , V. N. , Brindley , J. , Holden , A. V. and Tsyganov , M. A. 2004 . Pursuit-evasion predator-prey waves in two spatial dimensions . Chaos , 91 : 218102
- Britton , N. F. 1986 . Reaction-Diffusion Equations and Their Applications to Biology , London : Academic Press .
- Chakraborty , A. , Singh , M. , Lucy , D. and Ridland , P. 2007 . Predator-prey model with prey-taxis and diffusion . Math. Comput. Model. , 46 : 482 – 498 .
- Dolak , Y. and Hillen , T. 2003 . Cattaneo models for chemosensitive movement numerical solution and pattern formation . J. Math. Biol. , 46 : 153 – 170 .
- Garvey , J. E. , Stein , R. A. and Thomas , H. M. 1994 . Assessing how fish predation and interspecific prey competition influence a crayfish assemblage . Ecology , 75 ( 2 ) : 532 – 547 .
- Grunbaum , D. 1998 . Using spatially explicit models to characterize foraging performance in heterogeneous landscapes . Am. Nat. , 151 : 97 – 115 .
- Hamilton , W. D. 1971 . Geometry for the selfish herd . J. Theor. Biol. , 31 : 295 – 311 .
- Hill , A. M. and Lodge , D. M. 1999 . Replacement of resident crayfishes by an exotic crayfish: the roles of competition and predation . Ecol. Appl. , 9 ( 2 ) : 678 – 690 .
- Hillen , T. and Painter , K. 2001 . Global existence for a parabolic chemotaxis model with prevention of overcrowding . Adv. Appl. Math. , 26 ( 4 ) : 280 – 301 .
- Jorné , J. 1974 . The effects of ionic migration on oscillations and pattern formation in chemical systems . J. Theor. Biol. , 43 : 375 – 380 .
- Kareiva , P. and Odell , G. 1987 . Swarms of predators exhibit ‘preytaxis’ if individual predators use arearestricted search . Am. Nat. , 130 : 233 – 270 .
- Kolokolnikov , T. , Erneux , T. and Wei , J. 2006 . Mesa-type patterns in the one-dimensional brusselator and their stability . Phys. D , 214 : 63 – 77 .
- Kurganov , A. and Tadmor , E. 2000 . New high resolution central schemes for nonlinear conservation laws and convection-diffusion equations . J. Comput. Phys. , 160 : 214 – 282 .
- Lee , J. 2006 . Prey-taxis and its applications , Edmonton, , AB, Canada : University of Alberta . Ph.D. thesis
- Lee , J. , Hillen , T. and Lewis , M. A. 2007 . Prey-taxis equations: continuous travelling waves for pest control . in preparation
- LeVeque , R. J. 2002 . Finite Volume Methods for Hyperbolic Problems , Cambridge : Cambridge University Press .
- Lewis , M. A. 1994 . Spatial coupling of plant and herbivore dynamics: the contribution of herbivore dispersal to transient and persistent ‘waves’ of damage . Theor. Population Biol. , 45 : 277 – 312 .
- Losey , J. E. and Denno , R. F. 1998c . The escape response of pen aphids to foliar-foraging predators: factors affecting dropping behaviour . Eco. Entom. , 23 : 53 – 61 .
- Malchow , H. 1995 . Flow-and locomotion-induced pattern formation in nonlinear population dynamics . Ecol. Model. , 82 : 257 – 264 .
- Malchow , H. 2000 . Motional instability in prey-predator system . J. Theor. Biol. , 204 : 639 – 647 .
- Morton , T. L. , Haefner , J. , Nugala , V. , Decino , R. D. and Mendes , L. 1994 . The selfish herd revisited: do simple movement rules reduce relative predation risk . J. Theor. Biol. , 167 : 73 – 79 .
- Murray , J. D. 1989 . Mathematical Biology , New York : Springer Verlag .
- Okubo , A. and Levin , S. A. 2000 . Diffusion and Ecological Problems: New Perspectives , 2 , New York : Springer Verlag .
- Owen , M. R. and Lewis , M. A. 2001 . How predation can slow, stop or reverse a prey invasion . Bull. Math. Biol. , 63 : 655 – 684 .
- Rovinsky , A. B. and Menzinger , M. 1992 . Chemical instability induced by a differential flow . Phys. Rev. Lett. , 69 : 1193 – 1196 .
- Segel , L. A. 1980 . Mathematical Models in Molecular and Cellular Biology , Cambridge : Cambridge University Press .
- Segel , L. A. and Jackson , J. 1972 . Dissipative structure: an explanation and an ecological example . J. Theor. Biol. , 37 : 545 – 559 .
- Shi , W. and Zusman , D. R. 1993 . Fatal attraction . Nat. Sci. , 366 : 414 – 415 . Correspondence
- Strikwerda , J. C. 1989 . Finite Difference Schemes and Partial Differential Equations , CA : Wadsworth, Inc .
- Tsyganov , M. A. , Brindley , J. , Holden , A. V. and Biktashev , V. N. 2003 . Quasisoliton interaction of pursuitevasion waves in a predator-prey system . Phys. Rev. Lett. , 91 : 218102
- Turchin , P. 2003 . Complex Population Dynamics , Princeton and Oxford : Princeton University Press .
- Tyson , R. , Lubkin , S. R. and Murray , J. D. 1999 . Model and analysis of chemotactic bacterial patterns in a liquid medium . J. Math. Biol. , 38 : 359 – 375 .
- Tyson , R. , Stern , L. G. and LeVeque , R. J. 2000 . Fractional step methods applied to a chemotaxis model . J. Math. Biol. , 41 : 455 – 475 .
- Viscido , S. V. and Wethey , D. S. 2002 . Quantitative analysis of fiddler crab flock movement: evidence for ‘selfish herd’ behaviour . Anim. Behav. , 63 : 735 – 741 .
- Viscido , S. V. , Miller , M. and Wethey , D. S. 2001 . The response of a selfish herd to an attack from outside the group perimeter . J. Theor. Biol. , 208 : 315 – 328 .