Population dynamic consequences of fearful prey in a spatiotemporal predator-prey system

1.   Introduction

Ecologists have not paid much attention to the predator-induced fear effects until recently, as they have assumed that predator-induced stresses are acute and transitory ^[28] and positively related to direct effects. In the presence of a predator, prey shows various behavioural and morphological responses, including changes in habitat use ^[1,9,13], foraging ^[20,41], vigilance ^[5], agammaegation ^[4,2], sensitivity to initial conditions ^[40] and movement patterns ^[10,29], so that they become difficult to detect, encounter and capture. Due to these flexible responses of prey to a changing risk of predation, the predator may impact the prey demography more than just by direct killing ^[18,19,8].

Other demographic processes like reproduction, malnutrition and disease occur in large time scale. Traditionally, it has been mistaken that these processes are the result of food scarcity or parasitic infection ^[7,17]. Therefore, until some experimental studies ^[42,32], there was no reason for considering the predator-induced stresses in ecology.

Empirical studies have shown that indirect effects of predators can be comparable or even larger than the direct effects. In their experiment on spiny water fleas Bythotrephes longimanus with three species of zooplankton in Lake Michigan and Lake Erie, Pangle et al. ^[24] found that over six combinations of location and depth, indirect or risk effects on population growth rates were more than seven times larger than that of direct effects. In their seminal experiments, Peckarsky et al. ^[26], Schmitz et al. ^[29] and Nelson et al. ^[23] recognized the influence of fear by comparing the impacts of predatory and risk-only predators on prey demography ^[26], density ^[29] and population dynamics ^[23] respectively. They have demonstrated that risk effects can account for a large proportion of total impact of predators on prey ^[27]. In the presence of predators, snowshoe hare shifts to relatively safe but less profitable microhabitats which may hasten the decline and lengthen the recovery phase of the snowshoe hare cycle ^[14]. This habitat shift lowered the overall body condition of female hare, which ultimately cost their reproductive output. Zanette et al. ^[42] have done a field experiment on song sparrows during an entire breeding season by actively preventing the direct predation using the electric fences and seines and the predation risk was managed using the predator call playbacks. They measure the impact of fear on the demography and observed the 40% reduction in the offspring produced per year in the population of free-living song sparrows. Travers et al. ^[32] by their experiment demonstrated that female birds experiencing the frequent nest predation laid fewer eggs in their next nest. Thus, an emerging view is that mere presence of a predator may impact the prey demography more powerful than that of direct predation.

To incorporate these risk effects in predator-prey interactions, Wang et al. ^[37] have proposed a predator-prey model with the cost of fear. They have shown that the high levels of fear can stabilize the system and avoid the paradox of enrichment in ecosystems. Hopf-bifurcation with fear can be both supercritical and subcritical, whereas the model system without fear can have only supercritical Hopf-bifurcation. Wang and Zou ^[38] have discussed the fear effect in predator-prey interactions with adaptive avoidance of predators. Wang and Zou ^[39] investigated the reaction-diffusion predator-prey model with the cost of anti-predator behaviours. They have obtained the necessary and sufficient conditions of spatial pattern formation for various predator-prey functional responses. Complex dynamics of diffusive Leslie-Gower model with Holling type IV functional response and nonlinear harvesting have been studied by Upadhyay et al. ^[35]. They have discussed the Turing pattern formation, Hopf-bifurcation and stability of bifurcating periodic solutions. Upadhyay et al. ^[33] have analyzed the propagation of Turing patterns and periodic Travelling wave in spatial plankton model. Spatial Leslie-Gower type predator-prey model with time delay has been studied by Wang et al. ^[36]. Diffusion-driven instability, delay induced instability and various interesting patterns are also explored. Sun et al. ^[30] have investigated the spatial patterns through diffusion-driven instability in modified Leslie-Gower predator-prey model with Holling type Ⅱ functional response.

In the present work, we have proposed and analyzed a Leslie-Gower type predator-prey model incorporating the cost of fear. Study of several pairs of the interacting species, ranging from house sparrows and European sparrow hawk to mule deer and mountain lion ^[31] shows that the theoretical results obtained from Leslie-Gower predator-prey model system with Holling type Ⅱ functional response, closely approximate the practical reality ^[15] and appropriate for modeling the broad range of species. Spatial extension of the model has been considered here to investigate the influence of cost of fear on the spatial distribution of species, which was not considered in the previous studies made by Wang et al. ^[37] and Wang and Zou ^[38]. Extensive numerical simulations are performed to validate the theoretical findings and to explore different types of spatial patterns. Biological implication of results has also been discussed. In Section 2, we have presented the model system incorporating the cost of fear due to anti-predator behaviours. Boundedness, equilibria and stability analysis of the temporal model system are given in Section 3. In Section 4, spatially explicit model system has been presented. Stability analysis, Turing instability, the existence of Hopf-bifurcation, direction of Hopf-bifurcation and stability of bifurcating periodic solutions have also been discussed. Numerical simulation results are carried out in Section 5. In last Section, we have concluded the work by briefly summarizing the analytical and numerical findings.

2.   Model formulation

Consider a Leslie-Gower type predator-prey model with Holling type Ⅱ functional response, where population density of prey is denoted by $x(t)$ and predator, $y(t)$ at time $t$. The model system is given by the following system of ordinary differential equations

subjected to the positive initial conditions

The underlying assumptions of the model system (2.1) are as follows:

(ⅰ) Prey population is growing logistically in the absence of predation. The logistic growth of prey is composed of three parts: growth term $ax$, natural death term $-dx$ and density-dependent death term $-bx^2$ (due to intraspecific competition between the prey species).

(ⅱ) Predator $y$ is of generalist type whose most favorite food is prey $x$. Predator is consuming prey with Holling type Ⅱ functional response and growing logistically with prey dependent carrying capacity $\frac{x}{s}$.

(ⅲ) The factor $-\frac{sy}{x}$ tells that the rate of growth of predator population is limited and causes a decrease in the rate of increase of the predator population as $y$ increases.

All parameters $a, \; d, \; b, \;p, \;D, \;r, \;s\;$ are positive and their definitions are given in Table 1.

Field experiments ^[42] have suggested that fear may influence the prey reproduction. To incorporate the impact of fear effects, we modify the model system (2.1) by multiplying the growth term of prey by a factor $f(l, y)$, which accounts for the cost of anti-predator response due to fear. The parameter $l$ reflects the level of fear which derives the anti-predator behaviours in prey. It is assumed that scared prey forage less and left their newborn less protected, which in turn reduces the growth rate of prey. Therefore, it is reasonable to assume that

For convenience, we have taken the following particular form of term $f(l, y)$:

leading to the following system:

3.   Analysis of the temporal model system

3.1.   Boundedness and positivity

Boundedness and positivity of the solutions of the model system (2.4) have been presented in this subsection.

Lemma 3.1. The solutions of the model system (2.4) are positive and eventually bounded i.e. there exist $T \geq 0$ such that $x(t)< M_1$ and $y(t)< M_2,$ $\forall t \geq T.$

Proof. The phase portrait of the model system (2.4) is shown in Figure 1. The nullclines of the system are $C_1: y = \frac{-(p+Al)+\sqrt{(p+Al)^2+4pl(B-A)}}{2pl}$, $A = (d+bx)(D+x)$, $B = a(D+x)$ on which $\frac{dx}{dt} = 0;$ and $C_2:y = \frac{x}{s}$, on which $\frac{dy}{dt} = 0.$ $C_1$ and $C_2$ partition the first quadrant into four parts $D_1, \;D_2, \;D_3,$ and $D_4$. The intersection of $C_1$ and $C_2$ is the unique positive equilibrium point of the model system. Set $L_1 = \{(x, y)\; \vert\; x = M_1, \; 0\leq y \leq M_2$} and $L_2 = \{(x, y)\;|\; 0\leq x \leq M_1, \; y = M_2\}.$ The rectangular region consisting of $L_1, \; L_2, \; x-\text{axis}\; \text{and}\; y-\text{axis}\;$ as boundaries is denoted by $D_R$ ($= D_1 \cup D_2 \cup D_3 \cup D_4)$. Clearly, $D_R$ attracts all the trajectories starting in the first quardrant. Thus, solutions of the model system (2.4) are eventually bounded.

Next, we have discussed the positivity of solutions of the model system (2.4). For this, we have shown that any trajectory starting in first quadrant cannot reach the $y$-axis. To this end, we only need to prove that trajectory cannot arrive the $y$-axis in $D_2$. From a given point $(x_0, y_0)\in D_2$ the time of a trajectory running from $(x_0, y_0)$ to $C_1$ is denoted by $T_1$ and $T_2(N)$ is the time of trajectory running from $(x_0, y_0)$ to $\frac{x_0}{N}$, $N\in \mathbb{N}$, $N\geq 2$. We estimate the times $T_1$ and $T_2.$

Hence $T_1$ is finite. Now,

Since

therefore $\exists$ some $N_0 \in \mathbb{N}$ such that

Hence $T_2(N_0)>T_1.$ This shows that the time of trajectory running to the $y$-axis is far longer than that to $C_1$, i.e. trajectory runs into $D_3$ before it reaches the $y$-axis. From the properties of vector field shown in Figure 1, the trajectory cannot reach the $y$-axis in $D_3$. Therefore, any trajectory starting in the first quadrant cannot reach the $y$-axis.

From the above discussion, we know that there is no homoclinic or heteroclinic orbit in the domain $D_R$. Hence, it is proved.

3.2.   Equilibria analysis

Equilibrium points of the model system (2.4) are the intersection of zero growth rate isoclines (the prey zero growth rate isocline and the predator zero growth rate isocline) given by $xF_1(x, y) = 0$ and $yF_2(x, y) = 0$. It is easy to verify that the system (2.4) possesses two nonnegative equilibrium points:

(ⅰ) The predator free axial equilibrium point $E_1 = \left(\frac{a-d}{b}, 0 \right)$, provided $a>d$.

(ⅱ) The positive interior equilibrium point $E^*(x^*, y^*)$, where $(x^*, y^*)$ is the positive solution of the following system of equations

By solving the system (3.5), we obtain that $x^* = sy^*$ and $y^*$ is the root of the cubic equation

where

Now, the existence and uniqueness of the interior equilibrium point $E^*(x^*, y^*)$ can be proved with the help of simple algebraic manipulations. It is clear that the model system (2.4) has a unique positive equilibrium point if the equation (3.6) has unique positive root. The constant term $\xi_3 = \frac{-(a-d)D}{bs^2l}$ is negative provided $a>d$. Therefore, if $a>d$ is satisfied then by Descartes rule of sign cubic equation (3.6) has only one positive root $y^*$. Further, for this value of $y^*$, the corresponding value of $x^*$ is given by $x^* = sy^*$. From above discussion, it is clear that model system (2.4) posses only one positive interior equilibrium point $E^*(x^*, y^*)$.

Now, we have numerically illustrated the existence of the interior equilibrium point $E^*(x^*, y^*)$ (with the help of zero growth rate isoclines) and its dependence on the condition $a>d$ (cf. Figure 2). Parameter values are taken as $l = 0.01, \;b = 0.1, \;p = 1.1, \;D = 2.0, \;r = 0.16, \;s = 0.9.$

Remark 1. From the above discussion on equilibria analysis, it is clear that the model system (2.4) does not have any equilibrium point when $a<d$ (growth rate is less than the death rate of prey species). For $a>d$, both axial and the interior equilibrium point exist.

3.3.   Local stability analysis

Local stability of a model system indicates that the system is stable for short living perturbations.

Theorem 3.2. (i) The axial equilibrium point $E_1\left(\frac{a-d}{b}, 0 \right)$ is a saddle point having stable manifold along the $x$-axis and unstable manifold along the $y$-axis.

(ii) The positive interior equilibrium point $E^*(x^*, y^*)$ is locally asymptotically stable if the conditions $bx^*+r>\frac{px^*y^*}{(D+x^*)^2}$ and $\frac{px^*y^*}{(D+x^*)^2}<\frac{alx^*}{s(1+ly^*)^2}+bx^*$ hold.

Proof. For proof see the Appendix.

Example 3.1. For the set of biologically realistic parameter values (used in numerical simulations) $a = 1.2, \;l = 0.1, \;d = 0.1, \;b = 0.1, \;p = 1.1, \;D = 2, \;r = 0.16, \;s = 0.9$, the positive equilibrium point is obtained as $E^*(x^*, y^*) = (2.21257, 2.45841)$. We obtain that the conditions of Theorem 3.2, $bx^*+r = 0.381257>\frac{p x^* y^*}{(D+x^*)^2} = 0.337171$ and $\frac{p x^* y^*}{(D+x^*)^2} = 0.337171<\frac{alx^*}{s(1+ly^*)^2}+bx^* = 0.411326$ are satisfied. Also, $A_1 = 0.0440864>0, \;A_2 = 0.114577>0$, therefore, from Routh-Hurwitz criterion the positive equilibrium point $E^*(x^*, y^*)$ is locally asymptotically stable.

3.4.   Global stability of the interior equilibrium point

In this subsection, we have shown the global asymptotic stability of $E^*(x^*, y^*)$ by using the suitable Lyapunov function.

Theorem 3.3. Suppose the positive interior equilibrium point $E^*(x^*, y^*)$ is locally asymptotically stable then it is globally asymptotically stable.

Proof. Define a Lyapunov function

where $c>0$ has to be determined. Calculating the time derivative of $V$ along the solution of the model system (2.4) we have

As $m_1<x<M_1$ and $m_2<y<M_2$, thus we have

Therefore

where

From (3.9), it is obvious that $\frac{dV}{dt}<0.$ Consequently, except at $E^*$, $V$ is decreasing along all trajectories in the first quadrant. Thus, the model system (2.4) is globally asymptotically stable about $E^*(x^*, y^*)$.

4.   Spatially explicit model system

Species disperse in the space in search of food, shelter, mates and to evade predation. Therefore, the spatial element of ecological interactions is of main focus to understand how the ecological communities took shape. Ecological interactions involving the role of space are modeled with the help of reaction-diffusion systems. These systems of reaction-diffusion equations explain how the density of species distributed in the space affected by the local interactions of species and diffusion that causes the spread of species. It is assumed that the species are randomly spread in the domain $\Omega \subseteq \mathbb{R}^2$ so that Fick's law holds. Let $x(t, u, v)$ and $y(t, u, v)$ be the population densities of prey and predator respectively at time $t$ and position $(u, v)\in \Omega$. The spatiotemporal model corresponding to the temporal model (2.4) is given by the following reaction-diffusion system:

with the initial conditions

and the no-flux boundary conditions

where $d_1$ and $d_2$ are diffusion coefficients for prey and predator respectively, $\nu$ is the outward normal to $\partial\Omega$. The Laplace operator $\Delta = \frac{\partial^2 }{\partial u^2}$ in one-dimensional diffusion and $\Delta = \frac{\partial^2 }{\partial u^2}+\frac{\partial^2 }{\partial v^2}$ in two-dimensional diffusion.

4.1.   Stability analysis of the spatial model system

In order to discuss the linear stability analysis of the model system (4.1) about the spatially homogeneous steady state $E^*(x^*, y^*)$, we linearize this system using the transformations

where $\hat{x}(t, u, v)$ and $\hat{y}(t, u, v)$ are small time and space perturbations. Conventionally, $\hat{x}$ and $\hat{y}$ are taken as

where $0<\epsilon_1, \;\epsilon_2\ll 1$ and $\lambda_k$ is the wavelength. Also, $\mathtt{k} = (k_u, k_v)$ is the wave number vector and $k = |\mathtt{k}|$ is the wave number.

Substituting (4.4) and (4.5) in the spatial system (4.1) we have

where $a_{11} = -bx^*+\frac{px^*y^*}{(D+x^*)^2}$, $a_{12} = -\frac{px^*}{D+x^*}-\frac{alx^*}{(1+ly^*)^2}$, $a_{21} = \frac{r}{s}$, $a_{22} = -r$. The characteristic equation of the linearized system (4.6) is given by

where I is the identity matrix of order 2 and

In simplified form, the characteristic equation (4.7) can be rewritten as

where

Then, the solutions of the characteristic equation (4.8) yield the following dispersion relation

According to Routh-Hurwitz criterion for stability Re$(\lambda)<0$ if and only if

Theorem 4.1. Suppose that the interior equilibrium point $E^*(x^*, y^*)$ is locally asymptotically stable for the nonspatial model system (2.4). Then $E^*(x^*, y^*)$ is locally asymptotically stable for the spatial system (4.1) if the condition (4.10) holds.

Proof. The proof of the theorem directly follows from Routh-Hurwitz criterion, hence omitted.

4.2.   The existence of Hopf-bifurcation around the interior equilibrium point $E^*(x^*, y^*)$

Hopf-bifurcation is an instability induced due to the transformation of stability of focus ^[34]. In temporal systems, Hopf-bifurcation breaks the symmetry of the system and gives rise to the oscillations, which are periodic in time and uniform in space. Local Hopf-bifurcation for the spatial predator-prey model system has been studied by ^[35,43,34]. In this subsection, we have shown the existence of Hopf-bifurcation in the spatial predator-prey model system with fear effect. We have divided the whole analysis into the following points:

(ⅰ) Reduction of the model system to the origin $(0, 0)$

First, we translate the interior equilibrium point $E^*(x^*, y^*)$ to the trivial point using the transformation $\bar{x} = x-x^*, \;\; \bar{y} = y-y^*$. Then the model system (4.1) in the neighbourhood of origin can be rewritten as (bar signs are dropped out for simplicity of notations)

Using the Taylor series expansion about the $(x, y) = (0, 0)$, system (4.11) can be rewritten as

where

and

In the following, we use $h$ as the bifurcation parameter (in fact, $D$ is the bifurcation parameter, representing the extent of environmental protection to the prey). Therefore, in case of determining the stability of positive equilibrium solution $E^*(x^*, y^*)$ and existence of Hopf-bifurcation, variation in the value of $h$ plays an important role.

(ⅱ) Linearization of the reduced system about the origin $(0, 0)$

Now, we introduce new variables $X_1(t)$, $X_2(t)$ and $X(t)$ by

then the system (4.12) can be transformed into an abstract differential equation

where

and

Let

then

and

Now, linearizing the system (4.13) about the origin $(0, 0)$, one can obtain that

(ⅲ) The characteristic equation of the linearized system

The characteristic equation of the system (4.16) is given by

where

As we know that the stability of the trivial solution of (4.13) depends on the locations of roots of Eq. (4.17). When all roots of (4.17) have negative real parts, then the trivial solution of (4.13) is stable otherwise it is unstable.

Clearly, the eigenvalue problem

has eigenvalues $0 = \lambda_0<\lambda_1<\lambda_2<\cdots<\lambda_k<\cdots,$ and the corresponding eigenfunctions are

Let $\beta_{k}^{1} = \begin{pmatrix} \gamma_k \\ 0 \end{pmatrix}$ and $\beta_{k}^{2} = \begin{pmatrix} 0 \\ \gamma_k \end{pmatrix}$. Then $\mathcal{B}_k = \{\left(\beta_{k}^{1}, \beta_{k}^{2}\right)\}_{k = 0}^{\infty}$ construct a basis of the phase space of system (4.13). Therefore

can be decomposed as

Thus, the characteristic Eq. (4.17) is equivalent to

i.e.,

(ⅳ) Stability analysis and critical value of bifurcation parameter $h$

Let $J_k = K-\lambda_k D_M$ then

and

where

Thus, the characteristic equation (4.19) can be simply denoted as

According to (4.20) and det$K>0$, it is easy to know that under the condition (4.15), if $-(d_1+d_2)\lambda_k+h^*-h<0$ i.e. $h>h^*-(d_1+d_2)\lambda_k$ for every $k \in \mathbb{N}_{0}$, then the positive constant steady state solution $E^*(x^*, y^*)$ is stable; if there exist certain $k\in \mathbb{N}_{0}$ such that $h<h^*-(d_1+d_2)\lambda_k$ then the positive constant steady state solution $E^*(x^*, y^*)$ is unstable.

(ⅴ) Transversality condition

Suppose $i\omega$ is the pure imaginary root of equation (4.20), substituting $i\omega$ into the equation (4.20), we have

Separating real and imaginary parts of above equation, one can obtain that

Let

Then, the only $h$ at which the homogeneous Hopf-bifurcation occurs is $h = h_0$. Near $h_0$, substituting $\lambda = \theta_1+i\theta_2$ into (4.20) and separating real and imaginary parts, one can obtain that

Differentiating the Eq. (4.22) with respect to $h$, we get

Therefore, the transversality condition holds. According to the Hopf-bifurcation theorem for differential equations ^[6], we have the following results.

Theorem 4.2. Suppose $\left(bx^*+\frac{px^*}{s(D+x^*)}+\frac{alx^*}{s(1+ly^*)^2}\right)>\frac{p x^* y^*}{(D+x^*)^2}$.

Then

(i) If $h>h_0$, then the positive constant steady state solution $E^*$ of system (4.1) is locally asymptotically stable and unstable when $h<h_0$.

(ii) System (4.1) can undergoes a Hopf-bifurcation at the positive constant steady state solution $E^*(x^*, y^*)$ when $h = h_j$, $j \in \mathbb{N}_{0}.$

Example 4.1. To illustrate the results obtained in Theorem 4.2, we have numerically simulated the model system (4.1) using the pdepe Matlab tool. Parameter values are taken as $a = 1.2, \;l = 0.1, \;d = 0.1, \;b = 0.1, \; p = 1.1, \; r = 0.16, \; s = 0.90, \; d_1 = 0.001, \; d_2 = 0.001$. At $D = 2.0$, steady state solution is $E^*(x^*, y^*) = (2.21257, 2.45841 )$, $A_1 = 0.04407>0$, $A_2 = 0.11458>0$ and $bx^*+\frac{px^*}{s(D+x^*)}+\frac{a l x^*}{s(1+ly^*)^2} = 1.05327>\frac{px^*y^*}{(D+x^*)^2} = 0.33717$. We find $h = 0.15239>h_0 = 0.10830$. Thus, from Theorem 4.2, model system (4.1) is locally asymptotically stable. Again, at $D = 1.7$, $A_1 = -0.00103<0$, $A_2 = 0.11056>0$ and $bx^*+\frac{px^*}{s(D+x^*)}+\frac{a l x^*}{s(1+ly^*)^2} = 1.05828>\frac{px^*y^*}{(D+x^*)^2} = 0.36729$. Also, $h = 0.17807<h_0 = 0.17910$, therefore from Theorem 4.2, model system (4.1) is unstable about $E^*$ and a limit cycle appears in the small neighbourhood of the equilibrium point $E^*$. These facts are numerically illustrated in Figures 3 and 4.

4.3.   Direction and stability of the spatially homogeneous period orbits

In this subsection, we have derived the formulae for determining the direction and stability of bifurcated periodic solutions arising through Hopf-bifurcation using the center manifold theorem and normal form theory ^[12].

To determine the stability of bifurcated periodic solutions, we need to know the restriction of system to its center manifold at $h_0$. Denote by $L^*$ the conjugate operator

with domain

where the $H^2(\Omega)$ is the standard Sobolev space and

In fact, we choose

For all $\alpha \in D_{L^*}$, $\beta \in D_{L}$, it is easy to verify that

where $\langle \alpha, \beta\rangle = \int_{\Omega}\bar{\alpha}^{T}\beta dx$ denotes the inner product in $L^2(\Omega)\times L^2(\Omega)$.

Following Hassard et al. ^[12], we decomposed $X = X^{C}\oplus X^{S}$ with $X^{C}: = \left\lbrace zq+\bar{z} \bar{q}:z \in \mathbb{C} \right\rbrace,$ $X^{S}: = \left\lbrace W \in X: \langle q^*, W \rangle = 0 \right\rbrace.$ Then, for any $(x, y)^{T} \in X$, there exist a $z \in \mathbb{C}$ and $W = (W_1, W_2)^{T} \in X^{S}$ such that

which implies

Thus, the system in $(z, W)$ coordinates becomes

where $\tilde{f} = (f, g)^{T}$ and $f$, $g$ are defined by Eq. (4.12).

Now, the straightforward calculations show that

Notice that

On the center manifold, we have

and

This implies that

Therefore

where

According to Hassard et al. ^[12], we can obtain

With the help of above derived parameters, we can determine the direction of Hopf-bifurcation and stability of bifurcated periodic solutions. If $\mu_2>0(<0)$, then the Hopf-bifurcation is supercritical (subcritical); $\beta_2<0(>0)$, the bifurcated periodic solutions are stable (unstable).

4.4.   Turing instability

Turing instability breaks the spatial symmetry and leads to the formation of spatial patterns, which are stationary in time and oscillatory in space. Turing instability is also known as diffusion-driven instability as it arises due to the introduction of diffusion in the model system. For this instability to occurs in the spatial model system it is necessary that prey $x$ diffuses more slowly as compared to the predator $y$ so $d_1<<d_2$ is assumed.

Turing instability occurs about the homogeneous steady state $E^*(x^*, y^*)$, when the model system without diffusion (2.4) is stable and the model system with diffusion (4.1) is unstable for small perturbations around $E^*$

i.e.

Therefore, for diffusive instability to occurs in the spatial model system (4.1), it is necessary that the condition (4.10) fails to exist. However, it is evident that $\text{tr}(V_k)<0$, if the condition $A_1>0$ holds. Thus, we can not depend on the sign of $\text{tr}(V_k)$ for Turing instability. Therefore, diffusive instability occurs only when $\text{det}(V_k)<0$. Equivalently,

As $H(k^2)$ is a quadratic function of $k^2$, the minimum value of $H(k^2)$ is reached at $k^2 = k^2_{cr}$, where

Consequently, the condition for diffusion-driven instability is $H(k^2_{cr})<0$ i.e.

Theorem 4.3. The condition of Turing instability for system (4.1) is given by the conditions (4.28) and (4.29) leading to the condition

provided $E^*(x^*, y^*)$ is stable in the absence of diffusion.

Example 4.2. To illustrate the phenomenon of Turing instability numerically, we have taken the following set of parameter values

Plot of $H(k^2)$ as a function of $k$ is shown in Figure 5(a) for different values of level of fear ($l$). It is observed that interval of negativity decreases as the value of $l$ increases. Dispersion relation ($\text{Re}(\lambda$) vs. $k$) for different values of level of fear ($l$) is shown in the Figure 5(b). It is observed that interval of positivity for $\text{Re}(\lambda)$ is decreases as level of fear increases. Thus, occurrence of Turing instability is significantly affected by the value of level of fear ($l$).

Example 4.3. Space series of model system (4.1) at Turing parameter set (4.31) with $l = 0.1$ is given by Figure 6. It is observed that spatial distribution of species is irregular.

5.   Numerical simulations

5.1.   Numerical simulation results for the temporal model system

In this subsection, we have presented a series of numerical simulations to better understand the dynamics of the model system (2.4) and to support the analytical findings. Simulation experiments are carried out in Matlab (R2013a) using the fourth-order Runge-Kutta method. Base parameter values are chosen as

at which the model system (2.4) exhibits stable dynamics (cf. Figures 7(a) and 7(b)). This hypothetical set of parameter values is based on the values reported in Wang et al. ^[37] and Upadhyay et al. ^[35].

An emerging concept is that fear can significantly affect the prey dynamics. In Figures 8(a) and 8(b), we have plotted the population densities of prey and predator as a function of the level of fear ($l$). These figures indicate that prey and predator populations decrease with increase in the level of fear ($l$) and converge towards the minimum density as the level of fear ($l$) increases. Biologically, this implies that as the level of fear ($l$) increases, the anti-predator behaviors of prey cost the overall population size of prey. As the predator carrying capacity depends on prey density so corresponding predator population density also decreases. Furthermore, Figure 8 confirmed the theoretical argument that the higher values of prey anti-predator responses decrease the prey density as a result of the cost of fear.

In Figures 9(a) and 9(b), we have presented the time series plots of prey and predator with and without fear effects. The parameter set used here is given in Eq. (5.1) except $D = 1.8$. It is observed that, in the absence of fear the model system is showing the periodic behaviour and increasing the value of the level of fear ($l$) to $0.5$, dynamics becomes locally asymptotically stable. Biologically, this implies that risk effects have a stabilising impact on population dynamics of both prey and predator.

Bifurcation diagram taking the level of fear ($l$) as a bifurcation parameter is given in Figure 10(a). In this plot, we have taken $D = 1.0$ and other parameter values are same as given in Eq. (5.1). The supercritical Hopf-bifurcation is observed, as first Lyapunov coefficient is negative (cf. Figure 10(b)). Thus, there exists a stable limit cycle bifurcating from the equilibrium point $E^*(2.21257, 2.45841)$.

In Figures 11(a) and 11(b), we have presented the one parameter bifurcation diagrams of model system (2.4) taking predation rate ($p$) as the bifurcation parameter and other parameter values are same as given in Eq. (5.1). It has been observed that the model system is stable for $p<1.359287$ (both the eigenvalues of the system having the negative real parts). As we increase the values of $p$, a Hopf-bifurcation point $H(x = 1.811434, \;y = 2.012705, \;p = 1.359285)$ is obtained. At this point, first Lyapunov coefficient is $-2.689880e-02<0$ indicates a supercritical Hopf-bifurcation. Indeed, there are two eigenvalues of the equilibrium with $\text{Re}( \lambda)_{1, 2}\approx 0$ and $\text{Im}(\lambda)_{1, 2}\neq 0$ and first Lyapunov coefficient is negative. Thus, there exists a stable limit cycle at $p = 1.359285$ bifurcating from the equilibrium point (cf. Figure 11(c)). Further, increasing the value of $p$, we have observed that the model system is unstable for $1.359285<p<3.877143$ (both the eigenvalues of the system having the positive real parts) and we obtained another Hopf-bifurcation point as $H(x = 0.580594, \; y = 0.645104, \; p = 3.877143)$. At this point, first Lyapunov coefficient = $5.517181e-02>0$ indicates a subcritical Hopf-bifurcation. Indeed, there are two eigenvalues with $\text{Re}( \lambda)_{1, 2}\approx 0$ and $\text{Im}(\lambda)_{1, 2}\neq 0$ and first Lyapunov coefficient is positive. Thus, there exists an unstable limit cycle at $p = 3.877143$ bifurcating from the equilibrium point (cf. Figure 11(c)). In Figure 11(d), we have shown the period of cycle versus $p$ plot which is indicating the presence of two limit cycles with different periods for $p<4.475979$ near LPC (fold bifurcation of the cycle).

Figure 12(a) shows that prey are more willing to show anti-predator behaviors with increase in predation rate. Biologically, when predators more often attack prey (which cause more fear on prey species) and as a response to these fear effects, prey increases their vigilance and other anti-predator behaviors. Figure 12(b) shows that for higher values of environmental protection, prey species show less anti-predator responses, which support the theoretical argument that in relatively safer habitat and foraging sites prey species show fewer efforts towards fear.

5.2.   Simulation experiments for the spatial system

5.2.1.   1-D simulations

In this subsection, we have presented the simulation results for the model system (4.1) with one-dimensional diffusion. To exhibit the spatial dynamics, the system is solved numerically using the finite difference scheme. Forward difference scheme is used for the reaction part and central difference scheme is used for one-dimensional diffusion part. All the simulation experiments are performed in Matlab (R2013a). Simulations are carried out for different time and space steps until the obtained results become invariant. For simulation study, we assume spatial domain of the size 7000, space step as $\Delta h = 1$ and time step as $\Delta t = 0.001$. Parameter values and values of diffusion coefficient are chosen as

For realistic biological point of view, nonmonotonic form of initial conditions are assumed to describe the initial distributions of species $x$ and $y$ as

where $x^* = 2.21257$ and $y^* = 2.45841$ are prey and predator population densities at nontrivial steady state. $\epsilon = 10^{-8}$, $u_1 = 1200$ and $u_2 = 2800$ are parameters affecting the system dynamics.

Here, we have plotted the space vs. population density to understand the effect of movement of species in one-dimensional diffusion. Dynamics of prey and predator are observed at different time levels $t = 200, \; t = 500, \; t = 580$ and $t = 650,$ as shown in Figure 13. Initially, at $t = 200$ system is showing stable dynamics (cf. Figure 13(a)). As we increase time $t$ to 580, jagged patterns of small amplitude representing the irregular dynamics of the model system are obtained (cf. Figure 13(c)). These patterns grow slowly with time and their amplitude increases. At time $t = 650$, size of the whole domain is occupied by these irregular patterns (cf. Figure 13(d)). Chaotic attractor for the spatial system (4.1) with one-dimensional diffusion is presented in Figure 13(e). Thus, we can conclude that the spatial movement of species can introduce the irregular patterns in the otherwise stable system.

To investigate the spatiotemporal dynamics of the model system (4.1) with one-dimensional diffusion, we have presented Turing patterns as contour plots in Figure 14. In our simulations, initial conditions are specified as small random perturbations about the homogeneous steady state $E^*(x^*, y^*)$ with no-flux boundary conditions. Parameter values used in simulations are presented in Table 2. We have investigated the effect of the level of fear ($l$) on Turing patterns. At a relatively small value of the level of fear $l = 0.1$, prey and predator populations are distributed in the form of low and high-density patches as shown in Figures 14(a) and 14(b). Increasing the level of fear to $l = 0.5$, population density of prey and predator decreases, however nature of patterns remains unaltered (cf. Figures 14(c) and 14(d)). Again, increasing the level of fear to $l = 1.09$, prey and predator populations evenly distributed in the whole domain with low density as shown in Figures 14(e) and 14(f). Thus, biologically we can conclude that at a high value of fear, prey shows greater effort towards anti-predator responses which ultimately affect the predator-prey biomass and their distributions.

5.2.2.   2-D simulations

In this subsection, we have performed numerical simulations of the spatial model system (4.1) with two-dimensional diffusion. Reaction-diffusion system is solved by using the explicit Euler method for time integration and standard five-point approximation for the $2D$ Laplacian with no-flux boundary conditions. The model system is discretize in both space and time using the space step $\Delta h$ and time step as $\Delta t$. These step sizes are assumed sufficiently small to avoid numerical errors. This method leads to a sparse, banded linear algebraic system. The obtained linear system is solved using the GMRES algorithm ^[11].

Numerical simulations are carried out in Matlab (R2013a) with Turing parameter set given in Example 4.2. Size of the domain is assumed as $200\times 200$, time step $\Delta t = 0.001$, space step $\Delta h = 0.15$ and all numerical simulations employ Neumann boundary conditions. From Example 4.2 it follows that Turing instability occurs for those values of the fear level, which are less than $l = 1.09$. %and other parameters are given by (4.31).$l = 0.1, \; 0.5, \; 0.9$ holes-patches mixture pattern Now, we vary the value of the level of fear in Turing instability range and see the impact of it on the observed Turing patterns. At $l = 0.1$, random perturbation of the homogeneous steady state $E^*(2.21257, 2.45841)$ leads to the formation of holes pattern for both prey and predator (cf. Figures 15(a) and 15(b)). Increasing the value of fear level to $l = 0.5$, random perturbation of the homogeneous steady state $E^*(1.27340, 1.41489)$ leads to the formation of holes-small stripes mixture type of pattern for both prey and predator (cf. Figures 15(c) and 15(d)). Further, increasing the value of fear level to $l = 0.9$, random perturbation of $E^*(0.989014, 1.09890)$ leads to a significant decrease in the both prey and predator biomass (cf. Figures 15(e) and 15(f)).

Next, we have studied the impact of diffusion coefficient $d_2$ on Turing patterns with $l = 0.1$, $\Delta h = 0.05$, $E^* = (2.21257, 2.45841)$ and other parameters are taken from Turing space (see Example 4.2). When $d_2 = 0.1$, random perturbation of the homogeneous steady state $E^*$ leads to formation of holes pattern (cf. Figures 16(a) and 16(b)). Decreasing the value of $d_2$ to 0.09, observed Turing patterns remain unaltered (cf. Figures 16(c) and 16(d)). Further decreasing the value of $d_2$ to 0.05, prey and predator biomass decreases (cf. Figures 16(e) and 16(f)).

To obtain the nontrivial spatial patterns, initial distributions of species are considered as

where, $\epsilon_1 = 5\times10^{-4}$, $u_0 = v_0 = 0.1$, $S = 0.2.$

We have observed the evolution of Turing patterns for different values of level of fear at 1500 days. At $l = 0.1$, small random perturbation of the homogeneous steady state $E^*(2.21257, 2.45841)$ leads to the formation of holes-stripes mixture pattern for both prey and predator (cf. Figures 17(a) and 17(b)). Increasing the level of fear to $l = 0.5$, random perturbation of the homogeneous steady state $E^*(1.27340, 1.41489)$ tends to the formation of holes-small stripes mixed type pattern for both prey and predator (cf. Figures 17(c) and 17(d)).

In Figure 18, we have observed that holes-stripes mixture type pattern leads to the formation of holes pattern over the whole spatial domain by increasing the diffusion coefficient $d_1$ for prey.

6.   Discussions and conclusions

Belgrad and Griffen ^[3] investigated the predator-prey interactions in structural ecological communities mediated by prey personality and predator hunting mode. It has been found that the behaviour and survival of prey are affected by these factors which have a large potential to control trophic cascades and act as a mechanism for maintaining intraspecific trait variation. They may also guide predictions on the strength of predator-prey interactions as well as the response of ecosystems to such pervasive issues like species invasions, habitat destruction and overfishing. Various experimental studies on a wide range of species have suggested that the impact of fear on prey demography, density and dynamics can be large and sometimes can exceed the direct predation effects ^[7]. Zanette et al. ^[42] have reported the 40% reduction in the number of offspring of song sparrows due to fear of predator alone. Recently, some studies have been done to discover the impact of the cost of fear in predator-prey interactions involving the specialist predators ^[37,38,39].

In this work, we have tried to investigate the effect of fear on the more general scenario of predator-prey interactions involving generalist predators. We have formulated a Leslie-Gower type predator-prey model introducing the cost of fear on prey demography. The main focus of this study is to investigate the influence of anti-predator behaviours due to fear of predators in both space and time. Analytical results show that the positive coexistence equilibrium point $E^*(x^*, y^*)$ is locally as well as globally asymptotically stable, when certain conditions are satisfied. Turing instability, the existence of Hopf-bifurcation, direction and stability of bifurcating periodic solutions have been derived for the spatially explicit model system. Extensive numerical simulations are performed for both temporal and spatiotemporal model system. Variation in anti-predator responses with the change in other system parameters is demonstrated with the help of codimension two bifurcation diagrams. It has been observed that the increase in fear level can decrease both the prey and predator population densities (cf. Figure 8). Thus, a higher value of fear costs the overall population size. Bifurcation diagram with respect to the level of fear ($l$) is plotted (cf. Figure 10(a)). The supercritical Hopf-bifurcation phenomenon is observed. Bifurcation diagrams taking the predation rate $p$ as bifurcation parameter have also plotted, both the supercritical and subcritical Hopf-bifurcations have been observed (cf. Figure 11). It has been found that prey shows larger anti-predator behaviours for higher values of attack rate of predators (cf. Figure 12(a)). For environmentally protected shelters and foraging sites, prey are less willing to show anti-predator activities (cf. Figure 12(b)).

Diffusion has a destabilizing effect depending on the system parameters and nonlinearity ^[25]. By varying the level of fear, we have obtained the various Turing patterns for one and two-dimensional diffusive system. For instance, when $l = 0.1$, holes-stripes mixture type of patterns are observed (cf. Figures 17(a) and 17(b)). The significant decrease in prey and predator densities has been found for a further increase in the level of fear to $l = 0.5$ (cf. Figures 17(c) and 17(d)). We have also performed the numerical simulations with the fixed value of $d_2$ and by taking $d_1 = 0.001, \;0.003$ and $0.006$. Holes-stripes mixed type of patterns appeared which change to holes pattern by increasing the value of $d_1$ (cf. Figure 18).

In our study, we have considered the effect of fear in prey species and reported a significant change in population dynamics of both prey and predator with changing risk effects. Exposure to predation threat has been reported to incur significant impact on the fitness and phenotypes of animals across taxa ^[22,16]. McCauley et al. ^[21] showed the mere presence of a piscine predator, without access to prey, reduced dragonflies survivorship. Our study can add to increasing evidence of ecology of fear and have theoretical and experimental utility in the study and the regulation of crop pest in the ecoagricultural setting. To further our understanding of the nature of such a phenomenon we will investigate the differential suppression of prey reproduction and emergence of anti-predator defences in different model agroecosystems.

Appendix

The variational matrix of the model system (2.4) about an equilibrium point $E(x, y)$ is given by

let $V_{E_{1}}$ and $V_{E^*}$ denote the variational matrices of equilibrium points $E_1$ and $E^*$ respectively.

(ⅰ) At $E_1 = (\frac{a-d}{b}, 0)$

Eigenvalues about $E_1$ are $-a+d<0$ and $r>0$. As one eigenvalue is negative and other is positive, therefore axial equilibrium point $E_1$ is a saddle point having stable manifold along the $x$-axis and unstable manifold along the $y$-axis.

(ⅱ) At $E^*(x^*, y^*)$

The characteristic equation of above matrix is given by

let, $A_1 = -(a_{11}+a_{22})$, $A_2 = a_{11}a_{22}-a_{12}a_{21}.$

According to Routh-Hurwitz criteria $E^*(x^*, y^*)$ is locally asymptotically stable if

Thus, from Routh-Hurwitz criterion it is directly follows that $E^*(x^*, y^*)$ is locally asymptotically stable if the conditions

hold.

Acknowledgments

This work is supported by the SCIENCE & ENGINEERING RESEARCH BOARD(SERB), Department of science and Technology, Govt. of India, under File No. MTR/2017/000301 to the first and corresponding author (R. K. Upadhyay).

Conflict of interest

The authors declare there is no conflict of interest.