Drivers of pattern formation in a predator–prey model with defense in fearful prey

Abstract

Existence of predator is routinely used to induce fear and anxiety in prey which is well known for shaping entire ecosystem. Fear of predation restricts the development of prey and promotes inducible defense in prey communities for the survival. Motivated by this fact, we investigate the dynamics of a Leslie–Gower predator prey model with group defense in a fearful prey. We obtain conditions under which system possess unique global-in-time solutions and determine all the biological feasible states of the system. Local stability is analyzed by linearization technique and Lyapunov direct method has been applied for global stability analysis of steady states. We show the occurrence of Hopf bifurcation and its direction at the vicinity of coexisting equilibrium point for temporal model. We consider random movement in species and establish conditions for the stability of the system in the presence of diffusion. We derive conditions for existence of non-constant steady states and Turing instability at coexisting population state of diffusive system. Incorporating indirect prey taxis with the assumption that the predator moves toward the smell of prey rather than random movement gives rise to taxis-driven inhomogeneous Hopf bifurcation in predator–prey model. Numerical simulations are intended to demonstrate the role of biological as well as physical drivers on pattern formation that go beyond analytical conclusions.

1 Introduction

Fear of predation has immense impact on prey species. Predation directly increases the mortality of prey. Beyond the direct impact of mortality, fear of predation can elicit behavioral responses strong enough to affect the demography and life histories of entire prey populations [1]. Predation risk induces fear and can also influence the information-gathering strategies employed by prey. Fear results adaptation of a range of defensive behavior, which are aimed to escaping from the source of danger or motivational conflict [2]. Fear of predation indirectly affects on the growth of prey and also encourages prey to adopt anti-predation traits. Adoption of anti-predation defense is very common among prey. Some studies reported that social relationships play a key role in effectiveness and intensity of anti-predation defenses. Recently, Heathcote et al. [3] have reported that prey often joins larger groups in response to perceived predation risk or fear. For example, crested macaques (Macaca nigra) and dwarf mongooses (Helogale parvula) respond very quick and get associated with caller groups after getting predator alarm calls depending on their degree of affiliation with the caller [4]. Social prey usually forms herd to avoid or fight back with predator, whereas predator attacks on prey and kills them to take energy. But sometimes, killing prey is not easy to predators. Prey struggles for the survival and tries to avoid predation. Social prey prefers to live in group that helps them to avoid predation. Social prey makes group to fight back with predators that causes serious injury to the predators and reduces the predation rate. Some prey is even more dangerous to predators. Studies suggest that leopard attacks upon larger-sized prey, such as Cape buffalo (Syncerus caffer) and plains zebra (Equus quagga) that are too dangerous to hunt [5].

Recently, Brown and Laundre [6] have reported that anti-predator responses minimize risk of predation and help stabilize predator–prey dynamics in fear-driven systems. In mathematical ecology, the functional responses are used to describe the prey–predator relationship from cell to animals. Functional response describes the food availability per predator as a function of food density in per unit time. The functional response biologically considers handling time and attack rate of predators. Holling [7] proposed very realistic Holling disk equation for the predator–prey interaction of the form

$$\begin{aligned} f(u)=\frac{au}{1+a\tau u}. \end{aligned}$$

Holling considered handling time $\tau$ and encounter rate a in order to describe the relation between predator and their prey. Later on, some other functional responses are proposed of the form

• (I) Prey-dependent response:

  $f(x)=au$ (Holling I), $\displaystyle f(u)=\frac{au}{m+u}$ (Holling II),

  $\displaystyle f(u)= \frac{a u^2}{m^2+u^2}$ (Holling III).

• (II) Predator-dependent response:

  $\displaystyle f(u,v)=\frac{au}{\alpha u+\beta v+1}$ (Beddington–DeAngelis type),

  $\displaystyle f(u,v)=\frac{au}{(\alpha +u)(\beta +v)}$ (Crowley–Martin type)

• (III) Ratio-dependent: $\displaystyle f(u,v)=\frac{av}{u+v}$ (Michaelis–Menten type).

Holling I–III and all above-discussed responses are derived in the context of innocent prey and there is no component of defense or avoidance mechanism. Many biological studies reported that anti-predation usually increases the handling time because predator fights with individuals of large prey group. To model this phenomenon, we consider that predator handling time considers the time spent in fighting with every individual prey which is proportional to prey density. Let $t_f$ be the time spent in fight with a prey, then total spent time with fighting with each prey of group is $t_f u$. Thus, total time considered in case of group defense is $\tau =t_h+t_f u$. Thus, Holling disk equation can be rewritten as

$$\begin{aligned} f(u)=\frac{au}{1+at_h u+at_f u^2}, \end{aligned}$$

(1.1)

which considers the defense of prey against predators and known as Holling IV response [8]. A simplified version without handling time (i.e., $t_h=0$ in (1.1) $\displaystyle f(u)=\frac{au}{m+bu^2}$ of model (1.1) is proposed and applied in order to investigate group defense behavior in various predator–prey models [9,10,11]. On the basis of some field observations, it is reported that this non-monotonic function is best suitable functional response for anti-predation modeling [8, 12]. This functional response numerically supports the real-world biological data and also known as Monod–Haldane functional response [11].

Fig. 1

Holling functional responses

In Fig. 1, Holling I–III [7] are monotone function that describes the growth in food availability as food density increases. Holling IV is a non-monotonic function, i.e., the food availability starts decreasing after some critical values of food density. Biologically, this graph describes that predation is successful and predator can take more energy when group of prey is small. Moreover, rate of successful predation is decreased for large prey group size. Recently, Líznarová and Pekár [12] have shown that the non-monotonic consumption is associated with Holling IV response in dangerous prey such as spiders. From biological point of view, this phenomenon describe the group defense of social prey. An example of group defense in animals was firstly seen in muskoxen. Tener [13] observed that muskox make group to defend themselves from their predators wolf in Alaska. Another example of schooling-dependent defense is observed in small fish population of an aquatic system [14].

Recently, mathematical ecologists have proposed very specific mathematical models to study fear effect and/or group defense on predator–prey relationship with the consideration that fear has ability to reduce the growth of prey [15,16,17,18,19]. But according to our current information, only a ODE model [20] has been studied the simultaneous effects of fear and anti-predation. A very recent paper deals with a fractional order predator–prey model with fear and group defense [21]. However, there are several biological reports, suggesting that both mechanisms of fear and anti-predation occur during every encounter between prey and predator in every ecological complexity level. To our knowledge, this work is the first study that deals with diffusive predator–prey model with fear-induced group defense in prey. We report the emergence of coherent target wave, in general emerging due to Turing instability which is a novel outcome of this study. In this paper, we try to investigate the spatiotemporal dynamics of a diffusive predator–prey model. To study cascading effects of fear and group defense on moving predator–prey individuals, we propose the following predator–prey model

$$\begin{aligned} \left\{ \begin{aligned} \frac{\partial u }{\partial t}&=d_1 \varDelta u+ \frac{\alpha _1 u}{1+\kappa v}-\beta u^2-\frac{\gamma _1 uv}{bu^2+c}, ~~~~~~~~\mathbf{x } \in \varOmega ,\ t>0,\\ \frac{\partial v }{\partial t}&=d_2 \varDelta v+v\Big ( \alpha _2-\frac{\gamma _2 v}{m+u}\Big ) ~~~~~~~~~~~~~~~~~~~~~~~\mathbf{x } \in \varOmega ,\ t>0,\\ \frac{\partial u}{\partial \nu }&= \frac{\partial v}{\partial \nu }=0, ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\mathbf{x } \in \partial \varOmega ,\ t>0, \end{aligned} \right. \end{aligned}$$

(1.2)

with nonnegative initial distribution $u(\mathbf{x },0),\ v(\mathbf{x },0)>0$.

In model (1.2), the term $\displaystyle g(u,v,\kappa )= \frac{\alpha _1 u}{1+\kappa v}$ represents growth rate of prey as a function of predator density. Here $\displaystyle g(g,v,0)=\alpha _1 u,\ g(u,0,\kappa )=\alpha _1 u,\,\lim \nolimits _{\kappa \rightarrow \infty } g(u,v,\kappa )=0,\,\lim \nolimits _{v \rightarrow \infty } g(u,v,\kappa ){=}0, \frac{\partial g}{\partial \kappa }<0$ and $\displaystyle \frac{\partial g}{\partial v}<0$ suggest that development of prey individuals is suppose to decline when the fear increases or the number of predators increases. This fact is ubiquitous in ecology and is biological relevant to model fear effect in predator–prey systems. The parameter $\alpha _1$ measures the growth rate of prey and nonnegative parameter $\kappa$ represents the intensity of fear. The prey growth function $g(u,v,\kappa )$ without fear (i.e., $\kappa =0$) is $g(u,v)= \alpha _1 u$, which represents exponential growth in prey. We incorporate non-monotonic prey capture response to model group defense in predator–prey system with encounter rate $\gamma _1$. The parameter $\beta$ measures mortality rate due to intraspecific competition, b measures intensity of defense, c is saturation constant and $\alpha _2$ represents intrinsic growth rate of generalist predators. In order to deal with real-world situation, we assume that intraspecific competition among predators depends upon food availability and diffusion. The term $\displaystyle -\frac{\gamma _2 v^2}{m+u}$ represents intraspecific competition among predators depending upon availability of prey. This mathematical term of competition shows that food scarcity gives rise to high competition rate among predators. Here $\displaystyle \lim \nolimits _{u \rightarrow \infty } \frac{\alpha _2 v^2}{m+u}=0$ gives us an insight that there is no competition when food is abundant. This competition-induced mortality is common in predators which is studied by many authors and known as Leslie–Gower scheme [9, 22,23,24,25]. The nonnegative parameters $d_1$ and $d_2$ represent the rate of diffusivity of prey and predator individuals.

The study of proposed model is organized into 5 sections. Section 2 includes results on local- and global-in-time existence of solutions for the proposed model. We have discussed the properties of non-diffusive model in Sect. 3. Section 4 contains instability result of diffusive model. We have investigated the impact of chemical signaling on pattern formation in Sect. 5. Numerical simulation supporting theoretical results are discussed in Sect. 6. Finally we summarize the outcomes of our study in Sect. 7.

2 Existence of global positive solutions

Let us denote $x_+:x\vee 0, x_{-}:(-x)\vee 0$ for $x,y\in I\! \! R$ and $V_+=\{v\in V:v\ge 0\}$ which is positive cone of a partial order vector space $(V,\ge )$. Assume $Z=(u,v)\in X$ and its norm is defined as $\Vert Z\Vert =\Vert u\Vert _p+\Vert v\Vert _p$ where $p\in [1,\infty ),\ p<M/2$ and $X=\{L^p(\varOmega )\}^2$. The system (1.2) can be written as the following abstract equation:

$$\begin{aligned} \left\{ \begin{aligned}&Z_t=\nabla \cdot (A\nabla Z)+B(Z),~~~~~~~~~x\in \varOmega , t>0,\\&\frac{\partial Z}{\partial \nu }=0, ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ x\in \partial \varOmega , t>0,\\&Z(0;.)=(u_0(x),v_0(x)),~~~~~~~~~~~~~~ x\in \varOmega \end{aligned} \right. \end{aligned}$$

(2.1)

where

$$\begin{aligned} A&=\begin{pmatrix} d_1 &{} 0\\ 0 &{} d_2 \end{pmatrix},\ \ \nonumber \\ B(Z)&=\begin{pmatrix} u\Big (\displaystyle \frac{\alpha }{1+\kappa v}-\beta u-\frac{\gamma _1 uv}{bu^2+c}\Big )\\ v\Big (\displaystyle \alpha _2-\frac{\gamma _2}{m+u}\Big ) \end{pmatrix} \end{aligned}$$

(2.2)

and $Z=(u(t),v(t))^T$ with $\displaystyle \varLambda =\Big \{(u,v)\in \{C^2({\bar{\varOmega }}\}^2:\frac{\partial u}{\partial \nu }\Big |_{\partial \varOmega }=\frac{\partial v}{\partial \nu }\Big |_{\partial \varOmega }=0\Big \}$. Now we prove the results related to existence of local solutions for the system (2.1). Before discussing main results, we describe the following lemma:

Lemma 1

The system (2.1) with nonnegative initial condition (i.e., $Z_0\in X_+$) has a unique maximal local solution

$$\begin{aligned}&Z(t)\in ([0,t_{\max });X)\cap C^1((0,t_{\max });X)\nonumber \\&\times \cap C^1((0,t_{\max });\varLambda ), \end{aligned}$$

(2.3)

holds Duhamel formula for $t\in [0,t_{\max })$

$$\begin{aligned} Z(t)=e^{tA_p}Z_0+\int _{0}^t e^{(t-s)A_p}B[Z(s)]\mathrm{d}s, \end{aligned}$$

(2.4)

where operator $A_p$ represents a closure of $A\in X$. Moreover, if $t_{\max }<\infty$ then

$$\begin{aligned} \lim \sup _{t_{\max }} ||Z(t)||=\infty . \end{aligned}$$

Proof

The real part of eigenvalues of A(Z) is all positive, and hence, system (2.1) is uniformly parabolic system. Then, local existence, uniqueness of solutions and blowup follow from [26]. Here $A_p$ is closure of A in X, it generates a strong continuous and analytical semigroup $e^{tA_p}$. Furthermore, for $t>0$, we have

$$\begin{aligned} ||e^{tA_p}||\le e^{-t}||w||,\ \ \forall w\in X. \end{aligned}$$

Moreover, here we observe that $B:\varLambda \rightarrow X$ is locally Lipschitz on a bounded set. Thus, from Theorem 7.2.1 in [26], we can easily obtain desired result. $\square$

Lemma 2

For nonnegative initial data and for $p>N$, $(u_0,v_0)\in W^{1,p}(\varOmega )$, the solutions $(u(x,t),\ v(x,t))$ of the system (1.2) are also nonnegative.

Proof

The system (1.2) in positive cone can be written as

$$\begin{aligned} \left\{ \begin{aligned}&\frac{\partial {\bar{u}}}{\partial t}=d_1 \varDelta {\bar{u}}+ {\bar{u}}_+\Big ( \frac{\alpha _1 }{1+\kappa {\bar{v}}_+}-\beta {\bar{u}}_+-\frac{\gamma _1 {\bar{v}}_+}{b{\bar{u}}_+^2+c}\Big ), ~~~~~~~~~\mathbf{x } \in \varOmega ,\ t>0,\\&\frac{\partial {\bar{v}}}{\partial t}=d_2 \varDelta {\bar{v}}+{\bar{v}}_+\Big ( \alpha _2-\frac{\gamma _2 {\bar{v}}_+}{m+{\bar{u}}_+}\Big ) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\mathbf{x } \in \varOmega ,\ t>0, \end{aligned} \right. \end{aligned}$$

(2.5)

with initial and boundary conditions

$$\begin{aligned} \left\{ \begin{aligned}&\frac{\partial {\bar{u}}}{\partial \nu }= \frac{\partial {\bar{u}}}{\partial \nu }=0, ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\mathbf{x } \in \partial \varOmega ,\ t>0,\\&{\bar{u}}(x,0)=u_0,\ \ {\bar{v}}(x,0)=v_0 ~~~~~~~~~~~~~~~~~\mathbf{x }\in \varOmega . \end{aligned} \right. \end{aligned}$$

We multiply ${\bar{u}}_{-}$ and ${\bar{v}}_{-}$, respectively, in (2.5) and applying by part integration on domain $\varOmega$ we obtain

$$\begin{aligned} \left\{ \begin{aligned}&-\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}\int _{\varOmega } |\partial {\bar{u}}_{-}|^2 \mathrm{d}\mathbf{x }-d_1\int _{\varOmega } |\nabla \partial {\bar{u}}_{-}|^2 \mathrm{d}d\mathbf{x }=0,\\&-\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}\int _{\varOmega } |\partial {\bar{v}}_{-}|^2 \mathrm{d}\mathbf{x }-d_2\int _{\varOmega } |\nabla \partial {\bar{v}}_{-}|^2 \mathrm{d}\mathbf{x }=0. \end{aligned} \right. \end{aligned}$$

Therefore, for $t\in (0,t_{\max })$, we obtain

$$\begin{aligned} ||{\bar{u}}_{-}(.,t)||_2^2&+||{\bar{v}}_{-}(.,t)||_2^2\le ||{\bar{u}}_{-}(.,0)||_2^2\\&+||{\bar{v}}_{-}(.,0)||_2^2=0. \end{aligned}$$

Consequently, ${\bar{u}}(x,t)\ge 0,\ {\bar{v}}(x,t)\ge 0$. This proves the result. $\square$

It is known that $({\bar{u}}(x,t),{\bar{v}}(x,t))$ is the solution of the system (1.2) and Lemma (2) and we obtain $u(x,t)\ge 0,\ v(x,t)\ge 0,\ t\in (0,t_{\max })$. Now with the help of Lemmas (1) and (2), we are in the position to show global existence of positive solutions. To this end, it is sufficient to show that the solutions of (1.2) are uniformly bounded, i.e., dissipative. We follow lemma stated below:

Lemma 3

Suppose that p(x, t) is solution of the following equation

$$\begin{aligned} \left\{ \begin{aligned}&\frac{\partial p}{\partial t}=d\varDelta p+p(\alpha -p),~~~~~~~~~~~x\in \varOmega ,\ t>0\\&\frac{\partial p}{\partial \eta }=0 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~x\in \partial \varOmega ,\ t>0\\&p(x,0)=p_0(x)>0, ~~~~~~~~~~~~~~~~x\in \varOmega , \end{aligned} \right. \end{aligned}$$

Then $\lim _{t\rightarrow \inf } p(x,t)=\alpha$ for any $x\in \varOmega$.

Theorem 1

The solution (u(x, t), v(x, t)) of the system (1.2) satisfies the following inequalities

$$\begin{aligned} \lim \sup _{t\rightarrow \infty } \max _{x\in \varOmega } u(x,t)&\le \frac{\alpha _1}{\beta },\\ \lim \sup _{t\rightarrow \infty } \max _{x\in \varOmega } v(x,t)&\le \frac{\alpha _2 (\alpha _1 +m\beta )}{\gamma _2 m \beta }. \end{aligned}$$

Proof

From the equation of prey growth in (1.2),

$$\begin{aligned} \left\{ \begin{aligned}&\frac{\partial u}{\partial t}\le d_1\varDelta u+u(\alpha -\beta u),~~~~~~~~~~~~~~x\in \varOmega ,\ t>0\\&\frac{\partial u}{\partial \eta }=0 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~x\in \partial \varOmega ,\ t>0\\&p(x,0)=p_0(x)>0, ~~~~~~~~~~~~~~~~~~~~~~~~x\in \varOmega . \end{aligned} \right. \end{aligned}$$

From comparison inequality and Lemma (3), we can easily obtain first inequality. Let $\epsilon \in (0,1]$ is sufficiently small arbitrary constant then there exists a $T_1>0$ such that $\displaystyle u(x,t)\le \frac{\alpha _1}{\beta }+\epsilon$ for any $x\in \varOmega$ and $t>T_1$. Now for the second equation, we have

$$\begin{aligned} \left\{ \begin{aligned}&\frac{\partial v }{\partial t}\le d_2 \varDelta v\\&+ v\Big ( \alpha _2-\frac{\gamma _2 v}{\displaystyle m+\frac{\alpha _1}{\beta }+\epsilon }\Big ) ~~~~~~~~~~~~~~~~x\in \varOmega ,\ t>0\\&\frac{\partial v}{\partial \eta }=0 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~x\in \partial \varOmega ,\ t>0\\&p(x,0)=p_0(x)>0, ~~~~~~~~~~~~~~~~~~~~~~x\in \varOmega , \end{aligned} \right. \end{aligned}$$

From the comparison inequality and Lemma (3), for any $\epsilon \in (0,1]$, there exists at $T_2>T_1$ such that

$$\begin{aligned} v(x,t)\le \frac{\alpha _2 (\alpha _1 +m\beta )}{\gamma _2 m \beta }+2\epsilon . \end{aligned}$$

From the arbitraryness of $\epsilon$, we get

$$\begin{aligned} \lim \sup _{t\rightarrow \infty } \max _{x\in \varOmega } v(x,t)\le \frac{\alpha _2 (\alpha _1 +m\beta )}{\gamma _2 m \beta }, ~~~~~~ \forall x\in \varOmega . \end{aligned}$$

This completes the proof. $\square$

Lemmas 1 and 2 and Theorem 1 can be summarized as the following theorem

Theorem 2

The system (1.2) has an unique nonnegative and bounded solution $\big ((u(x,t),v(x,t)\big )$ such that

$$\begin{aligned} \big ((u(x,t),v(x,t)\big )\in C([0,\infty ),X)\cap C^1((0,\infty );X\\ \times \cap C((x,\infty ),\ \wedge (A)). \end{aligned}$$

3 Comments on ODE system

In this section, we assert stability and bifurcations in movement-free version of system (1.2). The model without movement factor can be written as set of nonlinear ordinary differential equations

$$\begin{aligned} \left\{ \begin{aligned} \frac{\mathrm{d} u }{\mathrm{d} t}&= \frac{\alpha _1 u}{1+\kappa v}-\beta u^2-\frac{\gamma _1 uv}{bu^2+c}, \\ \frac{\mathrm{d}v }{\mathrm{d} t}&=v\Big ( \alpha _2-\frac{\gamma _2 v}{m+u}\Big ),\\ \end{aligned} \right. \end{aligned}$$

(3.1)

with nonnegative initial condition $u(0)>0$ and $v(0)>0$. In the next subsection, we determine all the possible population states and conditions for their biological feasibility.

3.1 Stability of equilibrium points

In this subsection, we analyze the existence and stability of positive solution of ODE system (3.1). It is clear that the ODE system (3.1) has four equilibrium states. One can easily verify that $E_0=(0,\ 0)$ is trivial equilibrium of extinction. Predator extinction equilibrium $E_1=\displaystyle (\frac{\alpha _1}{\beta },\ 0)$ exists for positive values of $\alpha _1$ and $\beta$. Predator is generalist type, i.e., prey extinction equilibrium $E_2= \displaystyle (0,\ \frac{\alpha _2m}{\gamma _2})$ exists for positive parameters.

Now we discuss the existence of coexisting equilibrium state $E^\star =(u^\star ,v^\star )$. From the second equation of ODE system, we have

$$\begin{aligned} v^\star =\frac{\alpha _2(m+u^\star )}{\gamma _2}. \end{aligned}$$

(3.2)

Substituting $v^\star$ in first equation of ODE system, we get polynomial of the form

$$\begin{aligned} A_1 {u^\star }^4+A_2 {u^\star }^3 +A_3 {u^\star }^2+ A_4 {u^\star } +A_5=0, \end{aligned}$$

(3.3)

where

$$\begin{aligned}&A_1=\beta _1 \gamma _2 \alpha _2 b\kappa ,\ A_2=\beta \gamma _2 b(\alpha _2\kappa m+\gamma _2),\\&A_3=\alpha _2(\beta _1\gamma _2c\kappa +\alpha _2\gamma _1\kappa -b\gamma _2^2),\\&A_4=\gamma _2(\beta _1\alpha _2c\kappa m+\beta \gamma _1c+2\gamma _1\gamma _2\kappa m+\alpha _2 \gamma _2),\\&\ A_5=\alpha _1\kappa m^2+\alpha _2 \gamma _1 \gamma _2m-\alpha _2\gamma _2^2 c. \end{aligned}$$

In order to get at least a positive root we fix $A_5<0$ and choosing or $A_3<0$ or $A_3>0$ ensures at least a positive real root of $u^\star$. Hence, unique coexisting equilibrium $E^\star =(u^\star ,v^\star )$ exists under some certain restrictions. We numerically choose a set of parameters which gives a unique interior equilibrium $E^\star$. Now we examine the stability of each biological feasible equilibrium. We perturb the ODE system nearby its equilibrium point. For this, we calculate the Jacobian matrix at different equilibrium points.

• (i) Jacobian matrix calculated at $E_0=(0,0)$ is

  $$\begin{aligned} J_0=\begin{pmatrix} \alpha _1 &{} 0\\ 0 &{} \alpha _2 \end{pmatrix}. \end{aligned}$$

  Here both eigenvalues are positive. Hence, trivial equilibrium $E_0$ is always unstable.

• (ii) The Jacobian matrix calculated at $E_1=\displaystyle (\frac{\alpha _1}{\beta },\ 0)$

  $$\begin{aligned} J_1=\begin{pmatrix} -\alpha _1 &{}\displaystyle -\frac{\alpha _1\gamma _1\beta ^2+\alpha _1^2\kappa }{\beta (\alpha _1^2b+\beta ^2c)}\\ 0 &{} \alpha _2 \end{pmatrix}. \end{aligned}$$

  Jacobian matrix $J_1$ has two eigenvalues $-\alpha _1$ and $\alpha _2$ of the opposite signs. Hence, predator-free equilibrium $E_1$ is a saddle point.

• (iii) Jacobian matrix at $E_2= \displaystyle (0,\ \frac{\alpha _2m}{\gamma _2})$

  $$\begin{aligned} J_2=\begin{pmatrix} \displaystyle \frac{\alpha _1 \gamma _2}{\gamma _2+\kappa \alpha _2 m}-\frac{\gamma _1 \alpha _2 m}{c\gamma _2} &{} 0\\ \displaystyle \frac{\alpha _2^2 m^2}{\gamma _2} &{}\displaystyle -\alpha _2 \end{pmatrix}. \end{aligned}$$

  Here $\displaystyle -\alpha _2$ and $\displaystyle \frac{\alpha _1 \gamma _2}{\gamma _2+\kappa \alpha _2 m}-\frac{\gamma _1 \alpha _2 m}{c\gamma _2}$ are eigenvalues of matrix $J_2$. If $\displaystyle \frac{\alpha _1 \gamma _2}{\gamma _2+\kappa \alpha _2 m}<\frac{\gamma _1 \alpha _2 m}{c\gamma _2}$ holds, then $E_2$ is locally asymptotically stable and becomes unstable if it violated. Here one eigenvalue is negative. Thus, $E_2$ is a saddle point if $\displaystyle \frac{\alpha _1 \gamma _2}{\gamma _2+\kappa \alpha _2 m}>\frac{\gamma _1 \alpha _2 m}{c\gamma _2}$ holds. Now in particular $\displaystyle \frac{\alpha _1 \gamma _2}{\gamma _2+\kappa \alpha _2 m}=\frac{\gamma _1 \alpha _2 m}{c\gamma _2}$ holds; then, we have critical value of parameter $\alpha _1$ such that at

  $$\begin{aligned} \alpha _1=\alpha _1^\star =\frac{\gamma _1 \alpha _2 m(\gamma _2+\kappa \alpha _2 m)}{c\gamma _2^2} \end{aligned}$$

  the Jacobian matrix $J_2$ reduces to

  $$\begin{aligned} J_2\Big |_{\alpha _1=\alpha _1^\star }=\begin{pmatrix} \displaystyle 0 &{} 0\\ \displaystyle \frac{\alpha _2^2 m^2}{\gamma _2} &{}\displaystyle -\alpha _2 \end{pmatrix}. \end{aligned}$$

  (3.4)

  Clearly matrix (3.4) have a negative eigenvalue and a zero eigenvalue. If $V=(v_1,v_2)^T$ and $W=(w_1,w_2)^T$ denote eigenvectors corresponding to zero eigenvalue of $J_2$ & $J_2^T$ at ${\alpha _1=\alpha _1^\star }$, respectively, then we obtain $V=(1, \ \displaystyle \frac{\alpha _2m^2}{\gamma _2})^T$ and $W=(1,0)^T$. Now $\varDelta _1=W^T[f_{\alpha _1}(E_2,\alpha _1^\star )]=0$ suggests that saddle node bifurcation is not possible around $E_2$. Again in order to confirm transcritical bifurcation we calculate $\displaystyle \varDelta _2=W^T[Df_{\alpha _1}(E_2,\alpha _1^\star )V]=\frac{\gamma _2}{\gamma _2+\kappa \alpha _2 m}\ne 0$ where $\displaystyle Df_{\alpha _1}(E_2,\alpha _1^\star )V=\begin{pmatrix} \frac{\gamma _2}{\gamma _2+\kappa \alpha _2 m} &{} 0\\ 0 &{} 0 \end{pmatrix}.$ Similarly following [27], one can easily show $\displaystyle \varDelta _3=W^T[D^2f(E_2,\alpha _1^\star )(V,V)]\ne 0$ under the existence of the interior equilibrium point $E^\star$. Therefore, $E_2$ undergoes for a transcritical bifurcation at $\displaystyle \alpha _1^\star =\frac{\gamma _1 \alpha _2 m(\gamma _2+\kappa \alpha _2 m)}{c\gamma _2^2}$.

• (iv) Jacobian matrix at coexisting equilibrium point $E^\star$ is calculated as

  $$\begin{aligned} J=\begin{pmatrix} A_{11} &{}A_{12}\\ A_{21}&{}A_{22} \end{pmatrix}, \end{aligned}$$

  (3.5)

  where $\displaystyle A_{11}=\frac{2b\gamma _1{u^\star }^2v^\star }{(b{u^\star }^2+c)^2}-\beta u^\star$

  $$\begin{aligned}&A_{12}=-u^\star \Big (\frac{\gamma _1}{(b{u^\star }^2+c)}+\frac{\alpha _1\kappa }{(1+\kappa {v^\star }^2)}\Big ),\\&A_{21}=\frac{\gamma _2 {v^\star }^2}{(m+u^\star )^2},\ A_{22}=-\frac{\gamma _2 v^\star }{m+u^\star }. \end{aligned}$$

  The characteristic polynomial of $J_{E^\star }$ is

  $$\begin{aligned} |J-I \lambda |=\lambda ^2-\phi _1 \lambda +\phi _2=0, \end{aligned}$$

  (3.6)

  where $\phi _1=-(A_{11}+A_{22})$ and $\phi _2=A_{11}A_{22}-A_{12}A_{21}$. The Jacobian matrix J has negative eigenvalues if $\phi _1<0$ and $\phi _2>0$. It is easy to verify that $A_{22}<0$ and choosing $A_{11}<0$ or $\displaystyle \frac{2b\gamma _1{u^\star }^2v^\star }{(b{u^\star }^2+c)^2}<\beta u^\star +\frac{\gamma _2 v^\star }{m+u^\star }$ insure negative eigenvalues for the Jacobian matrix $J_{E^\star }$. Then $E^\star$ is locally asymptotically stable.

From the above analysis, we have the following theorem:

Theorem 3

The ODE system has the following stability properties

• (i) Trivial equilibrium $E_0$ is always unstable.

• (ii) Predator eradication state $E_1$ is saddle node.

• (iii) If (H1) $\displaystyle \frac{\alpha _1 \gamma _2}{\gamma _2+\kappa \alpha _2 m}<\frac{\gamma _1 \alpha _2 m}{c\gamma _2}$ holds, then prey eradication state $E_2$ is locally asymptomatic stable. $E_2$ becomes unstable if (H2) $\displaystyle \frac{\alpha _1 \gamma _2}{\gamma _2+\kappa \alpha _2 m}>\frac{\gamma _1 \alpha _2 m}{c\gamma _2}$ holds. $E_2$ undergoes for a transcritical bifurcation if $\displaystyle \frac{\alpha _1 \gamma _2}{\gamma _2+\kappa \alpha _2 m}=\frac{\gamma _1 \alpha _2 m}{c\gamma _2}$ holds.

• (iv) Either condition (H3) $2 b \gamma _1 u^\star v^\star <\beta (b{u^\star }^2+c)^2$ or (H4) $\displaystyle \frac{2b\gamma _1{u^\star }^2v^\star }{(b{u^\star }^2+c)^2}<\beta u^\star +\frac{\gamma _2 v^\star }{m+u^\star }$ holds, then coexisting equilibrium $E^\star$ is locally asymptotically stable.

Now we present the conditions for global stability $E^\star$ with the help of Lyapunov direct method. Global stability of solutions can be summarized by the following theorem:

Theorem 4

Suppose $E^\star$ is locally asymptotically stable and condition (A1) holds, then coexisting equilibrium $E^\star$ of ODE system (3.1) is globally asymptotically stable.

Proof

We apply Lyapunov direct method to examine the global stability of coexisting equilibrium $E^\star$ of ODE system (3.1). Let us consider the Lyapunov function of the form

$$\begin{aligned} L_1=u-u^\star -u^\star \log \big (\frac{u}{u^\star }\big )+\sigma \big [ v-v^\star -v^\star \log (\frac{v}{v^\star }\big )\big ] \end{aligned}$$

(3.7)

Derivative along time t gives

$$\begin{aligned} \frac{\mathrm{d}L_1}{\mathrm{d}t}=\frac{u-u^\star }{u}\frac{\mathrm{d}u}{\mathrm{d}t}+\sigma \frac{v-v^\star }{v}\frac{\mathrm{d}v}{\mathrm{d}t}. \end{aligned}$$

Here $\sigma$ is a positive parameter to be chosen suitably later. Some mathematical calculations give

$$\begin{aligned} \frac{\mathrm{d}L_1}{\mathrm{d}t}&=(u-u^\star )\Big [\alpha _1\Big (\frac{1}{1+kv}-\frac{1}{1+kv^\star }\Big )-\beta \big (u-u^\star \big )\\&-\gamma _1 \Big (\frac{v}{bu^2+c}-\frac{v^\star }{b{u^\star }^2+c}\Big )\Big ]\\&+\sigma (v-v^\star )\Big [ -\gamma _2 \Big (\frac{v}{m+u}-\frac{v^\star }{m+u^\star }\Big )\Big ]. \end{aligned}$$

Writing the above equation in the form of quadratics, we get

$$\begin{aligned} \frac{\mathrm{d}L_1}{\mathrm{d}t}&=-\Big [\beta -\frac{\gamma _1bv^\star (u+u^\star )}{(bu^2+c)(b{u^\star }^2+c)}\Big ](u-u^\star )^2\\&+\Big [ \frac{\sigma \gamma _2 v^\star }{(m+u^\star )}-\frac{\kappa \alpha _1}{(1+\kappa v)(1+\kappa v^\star )}-\frac{\gamma _1}{bu^2+c}\Big ]\\&\times (u-u^\star )(v-v^\star )-\frac{\sigma \alpha _2}{\delta _3} (v-v^\star )^2,\\ \frac{\mathrm{d}L_1}{\mathrm{d}t}&\le -\Big [\beta -\frac{\gamma _1bv^\star (u+u^\star )}{(b{u^\star }^2+c)}\Big ](u-u^\star )^2-\frac{\sigma \gamma _2}{m+u}\\&\times (v-v^\star )^2+\Big [ a_2\xi -\frac{\kappa \alpha _1 \gamma _2 m\beta }{\zeta _1}-\frac{\gamma _1\beta ^2}{b\alpha ^2+\beta ^2}\Big ]\\&\times (u-u^\star )(v-v^\star )\\, \end{aligned}$$

where $\displaystyle \zeta _1=(\kappa \alpha _1\alpha _2+m\beta (\gamma _2+\kappa \alpha _1))(1+\kappa v^\star )$. Now in particular, if we select $\displaystyle \xi =\frac{1}{a_2}\Big [\frac{\kappa \alpha _1 \gamma _2 m\beta }{\zeta _1}+\frac{\gamma _1\beta ^2}{b\alpha ^2+\beta ^2}\Big ]$, then we have

$$\begin{aligned} \frac{\mathrm{d}L_1}{\mathrm{d}t}&\le -\Big [\beta -\frac{\gamma _1bv^\star (u+u^\star )}{(b{u^\star }^2+c)}\Big ](u-u^\star )^2\\&\quad -\frac{\sigma \gamma _2}{m+u}(v-c^\star )^2,\\&=U^T D \ U, \end{aligned}$$

where $U^T=(u-u^\star ,\ v-v^\star )$, $\displaystyle D=\begin{pmatrix} \displaystyle \beta -\frac{\gamma _1bv^\star (u+u^\star )}{(b{u^\star }^2+c)} &{} 0\\ 0 &{}\displaystyle \frac{\sigma \gamma _2}{m+u} \end{pmatrix}.$. It is easy to verify that matrix D is a positive definite matrix if and only if

• (A1) $\displaystyle \beta >\frac{\gamma _1bv^\star (u+u^\star )}{(b{u^\star }^2+c)}$

holds. Hence, $\displaystyle \frac{\mathrm{d}L_1}{\mathrm{d}t}\le 0$ provided that $\displaystyle \beta >\frac{\gamma _1bv^\star (u+u^\star )}{(b{u^\star }^2+c)}$. Therefore, coexisting equilibrium $E^\star$ of ODE system (3.1) is globally asymptotically stable initiating in positive invariant set $\varLambda$. This proves our result. $\square$

3.2 Hopf bifurcation and its direction

In the previous section, we have investigated that the equilibrium points $e_0$ and $e_1$ are always unstable for positive parameters. Thus, there is no bifurcation at the equilibrium point $e_0$ and $e_1$. Some kind of bifurcation may take place at $e_2$ but we focus our attention on coexistence steady state $E^\star$. Now we examine the Hopf bifurcation near the coexisting equilibrium $E^\star$. In addition, we also discuss the direction of Hopf bifurcation at coexisting $E^\star$ when bifurcation parameter crosses its critical value $b=b^\star$.

Theorem 5

$E^\star$ undergoes for a Hopf bifurcation when parameter crosses its critical value. $b=b^\star$. Moreover, Hopf bifurcation is supercritical or subcritical depending upon $\rho <0$ or $\rho >0$, respectively.

Proof

The Jacobian matrix at $E^\star$ is presented in Eq. (3.6). When Jacobian matrix $J_{E^\star }$ has pure imaginary eigenvalue of the form $\lambda =\pm \omega$ and transversality condition holds, then Hopf bifurcation occur for $b=b^\star$. Obviously, it is sufficient to show that Hopf bifurcation occurs for some critical value $\alpha _{1}=\alpha _{1}^\star$, if the following conditions

• (i) $\displaystyle \phi _1\big |_{b^\star }=0$ and $\varPhi _2\big |_{b^\star }>0$,

• (ii) $\displaystyle \frac{\mathrm{d} \varPhi _1}{\mathrm{d} b}\big |_{b^\star }\ne 0$,

hold, where $\phi _1$ and $\phi _2$ are defined in Eq. (3.6).

Setting $\phi _1=0$ and performing some straightforward calculations, give the following threshold value

$$\begin{aligned} \alpha _{1}^\star =\frac{(\beta u^\star +\gamma _2v^\star )(b{u^\star }^2+c)}{2b(m+u^\star ){u^\star }^2v^\star } \end{aligned}$$

for which conditions (i) and (ii) hold. Now we focus on the direction of Hopf bifurcation and stability of periodic solutions staring nearby $E^\star$. In order to examine the stability of periodic solution, we calculate the Lyapunov coefficient at $\alpha _1=\alpha _1^\star$. We translate $E_0$ to the origin by transformation $u=u^\star -p$ and $v=v^\star -q$. Substituting these transformation in ODE model (3.1) and expanding with the help of Taylor series, we get

$$\begin{aligned} \left\{ \begin{aligned} \frac{\mathrm{d}p}{\mathrm{d}t}&=a_{10}p + a_{01}q + a_{20}p^2 + a_{11}pq + a_{02}q^2 \\&+a_{30}p^3+ a_{21}p^2q + a_{12}pq^2 + a_{03}q^3 + F(p,q),\\ \frac{\mathrm{d}q}{\mathrm{d}t}&=b_{10}p + b_{01}q + b_{20}p^2 + b_{11}pq + b_{02}q^2 \\&+ b_{30}p^3 + b_{21}p^2q + b_{12}pq^2 + b_{03}q^3 + G(p,q), \end{aligned} \right. \end{aligned}$$

(3.8)

where $\displaystyle F(p,q)=\sum \nolimits _{i+j=4}^{\infty } a_{ij} p^i q^j,\ \ \ \ \ F(p,q)=\sum \nolimits _{i+j=4}^{\infty } b_{ij} p^i q^j$,

and $a_{ij},\ b_{ij}$ are defined as

$$\begin{aligned} a_{ij}&=\Big [\frac{\partial ^{i+j}f_1(u,v)}{\partial u_i v_j}\Big ]_{\alpha _1=\alpha _1^\star },\ \ \\ b_{ij}&=\Big [\frac{\partial ^{i+j}f_2(u,v)}{\partial u_i v_j}\Big ]_{\alpha _1=\alpha _1^\star }. \end{aligned}$$

The first Lyapunov coefficient defined in [28] is given by

$$\begin{aligned} \rho&=-\frac{3\pi }{2a_{01}D^{3/2}}\Big \{ \big [a_{10}b_{10}(a_{11}^2 + a_{11}b_{02} + a_{02}b_{11})\\&\quad + a_{10}a_{01}(b_{11}^2 + a_{20}b_{11} + a_{11}b_{02})\\&\quad + b_{10}^2(a_{11}a_{02}+ 2a_{02}b_{02})-2a_{10}b_{10}(b_{02}^2-a_{20}a_{02})\\&\quad -2a_{10}a_{01}(a_{20}^2-b_{20}b_{02})-a_{01}^2(2a_{20}^2b_{20}+ b_{11}b_{20})\\&\quad + (a_{01}b_{10}-2a_{10}^2)(b_{11}b_{02}- a_{11}a_{20})\big ]-(a_{10}^2 \\&\quad + a_{01}b_{10})\big [3(b_{10}b_{03}-a_{01}a_{30})+2a_{10}(a_{21} + b_{12})\\&\quad +(b_{10}a_{12}-b_{01}b_{21})\big ]\Big ] \end{aligned}$$

For system (3.1), we have

$$\begin{aligned} a_{10}&=-\beta +\frac{2b\gamma _1u^\star v^\star }{(b{u^\star }^2+c)^2},\ \\ a_{01}&=-\Big [\frac{\gamma _1}{(b{u^\star }^2+c)}+\frac{\alpha _1 \kappa }{1+\kappa v^\star }\Big ],\\ a_{11}&=\frac{2b\gamma _1u^\star }{(b{u^\star }^2+c)^2}, \ a_{12}=0, \ a_{02}=\frac{2\alpha _1\kappa ^2}{(1+\kappa v^\star )^2}, \\ a_{20}&=\frac{2b\gamma _1v^\star (3b{u^\star }^2+c)}{b{u^\star }^2+c)^3},\ a_{21}=\frac{2b\gamma _1 (c-3b{u^\star }^2)}{(b{u^\star }^2+c)^3},\\ a_{03}&=-\frac{6\alpha _1\kappa ^3}{(1_+\kappa v^\star )^4},\ a_{30}=\frac{24b^2\gamma _1 u^\star v^\star (b{u^\star }^2-c) }{(b{u^\star }^2+c)^4},\\ b_{10}&=\frac{\gamma _2 v^\star }{(m+u^\star )^2},\ b_{01}=-\frac{\gamma _2 }{(m+u^\star )},\\ b_{11}&=-\frac{2\gamma _2 v^\star }{(m+u^\star )^2},b_{21}=\frac{2\gamma _2}{(m+u^\star )^3},\\ b_{12}&= b_{03}=b_{02}=0,\ b_{30}=\frac{6\gamma _2 v^\star }{(m+u^\star )^4} \end{aligned}$$

Let us consider $D=a_{10}b_{01}-a_{01}b_{10}>0$ and it is clear that $a_{10}+b_{01}=0$. From the above values, we have the following Lyapunov coefficient

$$\begin{aligned} \begin{aligned} \rho&=-\frac{3\pi }{2a_{01}D^{3/2}}\Big [a_{10}b_{10}(a_{11}^2+a_{02}a_{11})\\&+a_{10}a_{01}(b_{11}^2+a_{20}b_{11})\\&+b_{10}^2a_{11}a_{02}+2a_{10} a_{20} ( b_{10}a_{02}-a_{01}a_{20})\\&-a_{01}^2(2a_{20}b_{20}+b_{11}b_{20})-a_{11}a_{20}(a_{01}b_{10}-2a_{10}^2)\big ]\\&-(a_{10}^2+a_{01}b_{10})\big (2a_{10}a_{21}-a_{10}b_{21}-3a_{01}a_{30}\big ). \end{aligned} \end{aligned}$$

(3.9)

Fig. 2

Stable coexisting equilibrium $E^\star =(0.5685,\ 1.043)$ for $b=0.129<b^\star$ for the system (3.10)

The stability of periodic solution depends upon Lyapunov coefficient $\rho$. The limit cycle is stable if $\rho <0$ and unstable if $\rho >0$. Therefore, Hopf bifurcation is supercritical if $\rho <0$ and supercritical if $\rho >0$. $\square$

3.3 Numerical simulation of ODE system

Now we present some numerical results to illustrate the outcomes which are beyond the theoretic finding. Explicit formula and sign of Lyapunov coefficient is very tough to calculated. Thus, we calculate the sign of Lyapunov coefficient numerically in order to examine the stability and direction of Hopf bifurcation. To this end, we take the following example

$$\begin{aligned} \left\{ \begin{aligned} \frac{\mathrm{d} u }{\mathrm{d} t}&= u\Big (\frac{0.55}{1+\kappa v}-0.1 u-\frac{0.1 v}{b u^2+0.2}\Big ), \\ \frac{\mathrm{d}v }{\mathrm{d} t}&=v\Big ( 0.12-\frac{0.1 v}{0.3+u}\Big ),\\ \end{aligned} \right. \end{aligned}$$

(3.10)

with $u(0)=2.5,\ v(0)=3.0$. We present Figs. 2 and 4 as the evidence of Hopf bifurcation at $b=b^\star$. It is observed that system (3.10) has a stable positive equilibrium point $E^\star =(0.2299,\ 0.5160)$ for $b=0.162$ and $\kappa =0.1$. Theorem 4 is numerically illustrated in Fig. 3 which shows global stability of coexisting equilibrium E* According to Theorem 5, we show that the coexisting equilibrium becomes unstable for some $b>b^\star =0.165$. In Fig. 4, we observe that solutions of u(t) and v(t) oscillated around $E^\star$. This represents occurrence of limit cycle and destabilization of coexisting equilibria. Thus, Hopf bifurcation is verified for $b>b^\star$. In Fig. 5, we demonstrate the impact of fear on the dynamics of the system. It is observed that populations without fear and high defense efficiency oscillates around the coexistence population state. This oscillatory nature shows a stable limit cycle around $E^\star$ and demonstrates destabilization of coexistence population state in case of dangerous prey (see Fig. 5a). In Fig. 5b, it is observed that fear effect has an immense impact on the dynamics of the system. Fear in prey stabilizes the coexistence population state. However, fear has negative consequences on the evolution of prey but it has positive impact on the stability of the system. It is worth saying that the fear response is the opposite of the instability effect of prey group defense.

Fig. 3

Solutions converge to global attractor $E^\star =(0.5685,\ 1.043)$ for large perturbations

Fig. 4

Stable limit cycle surrounding the coexisting equilibrium $E^\star$ for $b=0.135>b^\star$ for the system (3.10)

Fig. 5

a Oscillatory populations around the coexistence steady state $E^\star$ without fear (i.e., $\kappa =0$), b populations converges to a constant steady state at $\kappa =0.13$ for the system (3.10)

Now in order to determine the stability and direction of Hopf bifurcation, we calculate coefficients $a_{ij}$ and $b_{ij}$ to obtain the sign of Lyapunov coefficient $\rho$. Numerical simulation gives the following at $b^\star =0.164$

$$\begin{aligned} a_{10}&=0.2278,\ a_{01}=-1.0112,\ a_{11}=0.6353,\ \\ a_{20}&=0.9735,\\ a_{02}&=0.0328,\ a_{21}=1.8866,\\ a_{12}&=0.0,\ a_{03}=-0.0178,\ a_{30}=-4.9669,\\ b_{10}&=0.2792, b_{01}=-0.2326,\ b_{11}=0.5411,\\ b_{20}&=-1.2989,\ b_{21}=-2.5172,\ b_{30}=9.0642. \end{aligned}$$

We obtain $\rho =-1.5723$ by substituting the above values in Eq. (3.9). Obviously, the Lyapunov coefficient is negative. Thus, following Theorem 5, we conclude that Hopf bifurcation is supercritical, i.e., limit cycle is stable. The conclusion follows results presented in Fig. 4.

4 Analysis of PDE system

In this section, we analyze the local and global stability of coexistence steady states in the presence of diffusion. The motive of this section is to investigate the impact of diffusion on the constant steady state. Apart from this, we focus on the situations which may give rise to instability. We discuss existence of Hopf and diffusion-driven instability for the diffusive model (1.2). We also present conditions for existence of non-constant steady state solution and discuss taxis-driven instability for diffusive system (1.2).

4.1 Stability of steady state

In this section, we derive the Turing bifurcation conditions around the homogeneous equilibrium $E^\star$ of the diffusive predator–prey model (1.2). For this purpose, we linearize the system by small perturbation nearby its steady state. We take the following form of perturbation near the equilibrium $E^\star$

$$\begin{aligned} \begin{pmatrix} u\\ v \end{pmatrix}= \begin{pmatrix} c_1^k\\ c_2^k \end{pmatrix}e^{\lambda _kt+ikr}, \end{aligned}$$

where k is wave number, $\lambda _k$ represents change in perturbation in time t and $c_2^k$ and $c_2^k$ are the amplitudes. Substituting the above transformations in Eq. (4.2), we get the following linearized form

$$\begin{aligned} \begin{aligned} \displaystyle \frac{\partial u}{\partial t}&=a_{11}u+a_{12}v-k^2 d_1 u,\\ \displaystyle \frac{\partial v}{\partial t}&=a_{21}u+a_{22}v-k^2 d_2 v. \end{aligned} \end{aligned}$$

(4.1)

The standard linear form of the PDE system (1.2) can be written as

$$\begin{aligned} \frac{\partial }{\partial t} \begin{pmatrix} u\\ v \end{pmatrix}= L\begin{pmatrix} u\\ v \end{pmatrix}:= D \begin{pmatrix} u_{xx}\\ v_{xx} \end{pmatrix}+\mathbf{J } \begin{pmatrix} u\\ v \end{pmatrix}. \end{aligned}$$

(4.2)

Here $D=\text {diag}(d_1,d_2)$ and L represents a linear operator defined in $H^2[\varOmega ]$ standard Sobolev space with domain $D_L=X_c:=X\oplus iX=\{x_1+ix_2: x_1,x_2\in X\}$ such that

$$\begin{aligned} X:=\Big \{(u,v)\in H^2[\varOmega ]\times H^2[\varOmega ]\Big | \begin{array}{l} \frac{\partial u}{\partial n}\big |_0=\frac{\partial u}{\partial n}\big |_l=0 \\ \frac{\partial v}{\partial n}\big |_0=\frac{\partial v}{\partial n}\big |_l=0 \end{array} \Big \}. \end{aligned}$$

J represents Jacobian matrix at steady state $E^\star$ and defined as

$$\begin{aligned} \mathbf{J }=\begin{pmatrix} \theta _0 &{} \sigma \\ s &{} -\theta _1 \end{pmatrix}, \end{aligned}$$

(4.3)

where $\theta _0= A_{11},\ \sigma =A_{12},\ s=A_{21}$ and $\theta = -A_{22}$ are defined in (3.5). From the above Jacobian matrix, we get

$$\begin{aligned} J_{s}=\begin{vmatrix} \theta _0-d_1k^2-\lambda&\sigma \\ s&-\theta _1-d_2k^2-\lambda \end{vmatrix}=0. \end{aligned}$$

The characteristic polynomial of the above Jacobian matrix is

$$\begin{aligned} {\lambda _k}^2-\varPhi _1(k^2){\lambda _k}+\varPhi _2 (k^2)=0, \end{aligned}$$

(4.4)

where

$$\begin{aligned} \varPhi _1(k^2)&=-(d_1+d_2)k^2+(\theta _0-\theta _1),\\ \varPhi _2(k^2)&=d_1d_2k^4+(d_1\theta _1-d_2\theta _0)k^2-(\theta _0\theta _1+\sigma s). \end{aligned}$$

The roots of (4.4) are

$$\begin{aligned} \lambda _k^{1,2}=\frac{\varPhi _1\pm \sqrt{\varPhi _1^2-4\varPhi _2}}{2}. \end{aligned}$$

(4.5)

If ($H_3$) holds then $\theta _1<0$. One can easily verify that $\sigma <0,\ s>0$ and $\theta _1>0$. Then $\varPhi _1(k^2)<0$ and $\varPhi _2(k^2)>0$, we conclude that the characteristic polynomial (4.4) have roots with negative real part for all $k\ge 0$. Hence, if $\frac{d_1}{d_2}>\frac{\theta _0}{\theta _1}$ i.e predator is much faster than prey and conditions $(H_3)$ and ($H_4$) holds then $\varPhi _2(k^2)$ is positive. Therefore, the coexistence steady state $E^\star$ with diffusion is locally stable. The above analysis can be summarized by the following theorem:

Theorem 6

Assume that $E^\star$ is locally asymptotically stable without diffusion and $\displaystyle \frac{d_1}{d_2}>\frac{\theta _0}{\theta _1}$ holds in the presence of diffusion, then the coexistence steady state $E^\star$ with diffusion is also locally asymptotically stable for all $k\in {\mathbb {N}}_{+}$.

Thanks to Theorem (3.3) in [29], we are in a position to comment on the global stability of the system. For the similar Lyapunov function defined in (3.7) and the following analysis performed [29], we have the following conclusion:

Theorem 7

If coexistence steady state $E^\star$ of ODE version of model (1.2) is globally asymptotically stable, then $E^\star$ of diffusion-induced model (3.1) is also globally asymptotically stable.

4.2 Nonexistence of positive non-constant steady states

In this subsection, we present the conditions for nonexistence of non-constant steady state solution. To this end, we define

$$\begin{aligned} f(u,v)&=\frac{\alpha _1 u}{1+\kappa v}-\beta u^2-\frac{\gamma _1 uv}{bu^2+c} \ g(u,v) \nonumber \\&=v\Big ( \alpha _2-\frac{\gamma _2 v}{m+u}\Big ), \end{aligned}$$

(4.6)

Theorem 8

Suppose $d_2 \mu _1 >\alpha _2$. Then there exists a positive constant $d=d(\varLambda ,d_2^\star )$ such that (1.2) has no non-constant positive solution provided that $d_2\ge d_2^\star$.

Proof

Let us assume that (u, v) is a positive solution of (1.2). Let ${\bar{\varphi }}= \frac{1}{|\varOmega |}\int _{\varOmega } \varphi \mathrm{d}x$ for $\varphi \in L^1(\varOmega )$. We multiply $u-{\bar{u}}$ to the first equation of (1.2) and then integrating on $\varOmega$

$$\begin{aligned}&d_1 \int _{\varOmega } |\nabla (u-{\bar{u}})|^2 \mathrm{d}x\\&\quad =\int _{\varOmega } f(u,v) (u-{\bar{u}}) \mathrm{d}x,\\&=\int _{\varOmega } \{f(u,v)-f({\bar{u}},{\bar{u}})\} (u-{\bar{u}}) \mathrm{d}x,\\&=\int _{\varOmega } \{f_u(\xi ,\mu ) (u-{\bar{u}})^2+f_v(\xi ,\mu )(u-{\bar{u}}) (v-{\bar{v}})\} \mathrm{d}x,\\&\le \int _{\varOmega } \{c_1 (u-{\bar{u}})^2+c_2 |(u-{\bar{u}})| |(v-{\bar{v}})| \mathrm{d}x. \end{aligned}$$

Here $c_1$, $\epsilon$ are positive constants and $\xi \in (u,{\bar{u}})$ and $\mu \in (v,{\bar{v}})$. Similarly we have

$$\begin{aligned}&d_2 \int _{\varOmega } |\nabla (v-{\bar{v}})|^2 \mathrm{d}x\\&=\int _{\varOmega } g(u,v) (v-{\bar{v}}) \mathrm{d}x,\\&=\int _{\varOmega } \{g(u,v)-g({\bar{u}},{\bar{u}})\} (v-{\bar{v}}) \mathrm{d}x,\\&=\int _{\varOmega } \{g_u(\xi ,\mu ) |(u-{\bar{u}})| |(v-{\bar{v}})|\\&\quad +\Big (\alpha _2-\frac{\gamma _2}{m+u}\Big )(v-{\bar{v}})^2\} \mathrm{d}x,\\&\le \int _{\varOmega } \{c_3 |(u-{\bar{u}})| |(v-{\bar{v}})|+\alpha _2 (v-{\bar{v}})^2\}\mathrm{d}x. \end{aligned}$$

Adding the above quantities, we get

$$\begin{aligned}&\int _{\varOmega } (d_1|\nabla (u-{\bar{u}})|^2 + d_2 |\nabla (v-{\bar{v}})|^2) \mathrm{d}x \\&\le \int _{\varOmega } \{c_1 (u-{\bar{u}})^2+2 \nu |(u-{\bar{u}})| |(v-{\bar{v}})|\\&\quad +\alpha _2 (v-{\bar{v}})^2 \mathrm{d}x\\&\le \int _{\varOmega } \{(u-{\bar{u}})^2\Big (c_1+\frac{\nu }{\epsilon }\Big )+ (v-{\bar{v}})^2 (\alpha _2+\nu \epsilon )\mathrm{d}x\}, \end{aligned}$$

for positive constant nu and some arbitrary constant $\epsilon >0$. Applying Poincáre inequality, we obtain

$$\begin{aligned}&\mu _1 \int _{\varOmega } (d_1|\nabla (u-{\bar{u}})|^2+d_2 |\nabla (v-{\bar{v}})|^2 )\mathrm{d}x \\&\le \int _{\varOmega } \{(u-{\bar{u}})^2\Big (c_1+\frac{\nu }{\epsilon }\Big )+ (v-{\bar{v}})^2 (\alpha _2+\nu \epsilon ) \}\mathrm{d}x. \end{aligned}$$

Since $\mu _1 d_2>\alpha _2$, we may choose sufficiently small $\epsilon _1$ such that $d_2\mu _1>\alpha _2+\nu \epsilon _1$, Consequently, taking $d:=(1/\mu _1)(1+\nu /\epsilon _1)$ implies that $u={\bar{u}}$ and $v={\bar{v}}$. This completes the proof. $\square$

4.3 Existence of positive non-constant steady states

In this subsection, we discuss the existence of non-constant steady state solutions depending upon diffusivity of individuals. For the simplicity, we denote $\mathbf{u }=(u,v)^T,\ \mathbf{F }=(f,g)^T,\ \mathbf{F }_u=J$ where J is defined in (4.3), and f and g are defined in (4.6). The system (1.2) can be written as

$$\begin{aligned} \left\{ \begin{aligned} -\varDelta \mathbf{u }&=D^{-1} \mathbf{F}(u) , ~~~~~~~~~~ x\in \varOmega ,\\ \frac{\partial \mathbf{u }}{\partial n}&=0 ~~~~~~~~~~~~~~~~~~~~~~ x\in \partial \varOmega , \end{aligned} \right. \end{aligned}$$

(4.7)

where $D=\text {diag}(d_1,d_2)$. Therefore, $u*$ solve (1.2) if and only if

$$\begin{aligned} {\bar{f}}(d_1,d_2,u)&=\mathbf{u }-(\mathbf{I }-\varDelta )^{-1}\{D^{-1}{} \mathbf{F}(u) +\mathbf{u }\}\nonumber \\&=0 \ \text {on} \ \mathbf{X }, \end{aligned}$$

(4.8)

where $(\mathbf{I }-\varDelta )^{-1}$ is inverse of $(\mathbf{I }-\varDelta )$ with homogeneous no-flux boundary condition. Straightforward calculations give

$$\begin{aligned} D_\mathbf{u }{\bar{f}}(d_1,d_2,u)=\mathbf{I }-(\mathbf{I }-\varDelta )^{-1}(D^{-1}J_{E^\star }+\mathbf{I }). \end{aligned}$$

For $\forall \mathbf{X }_i$ if and only of $\lambda$ is eigenvalue of the following matrix

$$\begin{aligned} M_i:=\mu _i I-D^{-1} J=\begin{pmatrix} \mu _i-d_1^{-1}\theta _0 &{} -d_1^{-1}\sigma \\ -d_2^{-1}s &{} \mu _i+d_2^{-1}\theta _1 \end{pmatrix}, \end{aligned}$$

where $\theta _0=A_{11},\ \sigma =A_{12},\ s=A_{21},\ \theta _1=-A_{22}$ defined in (3.5). It is easy to calculate $\det M_i=d_1^{-1} d_2^{-1}[d_1d_2\mu _i^2+(d_1\theta _1-d_2\theta _0)\mu _i-(s\sigma +\theta _0\theta _1)]$ and $\text {tr}\ M_i= 2\mu _i+d_2^{-1}\theta _1-d_1^{-1}\theta _0$. Let us define

$$\begin{aligned} {\bar{g}}(d_1,d_2,\mu )&=d_1d_2\mu ^2+(d_1\theta _1-d_2\theta _0)\mu \\&\quad -(s\sigma +\theta _0\theta _1), \end{aligned}$$

then ${\bar{g}}(d_1,d_2,\mu )=d_1d_2\det M_i$. If

$$\begin{aligned} (d_1\theta _1-d_2\theta _0)>-4d_1d_2(s\sigma +\theta _0\theta _1), \end{aligned}$$

(4.9)

then ${\bar{g}}(d_1,d_2,\lambda )=0$ has the following real eigenvalues

$$\begin{aligned}&\mu _{+}(d_1.d_2)=\frac{(d_2\theta _0-d_1\theta _1)+\sqrt{(d_1\theta _1-d_2\theta _0)^2+4d_1d_2(s\sigma +\theta _0\theta _1)}}{2d_1d_2},\\&\mu _{-}(d_1.d_2)=\frac{(d_2\theta _0-d_1\theta _1)-\sqrt{(d_1\theta _1-d_2\theta _0)^2+4d_1d_2(s\sigma +\theta _0\theta _1)}}{2d_1d_2}. \end{aligned}$$

Let us define

$$\begin{aligned} \varPi&=\varPi (d_1,d_2)=\{\mu :\mu \ge 0,\mu _{-}(d_1,d_2)<\mu \\&\qquad \qquad \qquad \qquad <\mu _{+}(d_1,d_2)\},\\ S_{p}&=\{\mu _0,\mu _1,\mu _2....\}, \end{aligned}$$

where $S_p$ is positive spectrum of $\varDelta$ on $\varOmega$ with homogeneous NBC. We apply the following lemma to calculate index of ${\bar{f}}(d_1,d_2;\cdot )$ at $E^\star$

Lemma 4

Suppose ${\bar{g}}(d_1,d_2;\mu _i)\ne 0,\ \forall \mu _i\in S_p$. Then $\text {index}({\bar{f}}(d_1,d_2;\cdot )=(-1)^\sigma$, where

$$\begin{aligned} \sigma = {\left\{ \begin{array}{ll} \sum _{\mu _i\in \varPi \cap S_p} m(\mu _i) &{}\quad \text {if} \quad \varPi \cap S_p\ne \emptyset \\ 0 &{}\quad \text {if}\quad \varPi \cap S_p=\emptyset . \end{array}\right. } \end{aligned}$$

In particular, if ${\bar{g}}(d_1,d_2;\mu _i)> 0$ for all $\mu _i\ge 0$, then $\sigma =0$.

Theorem 9

Suppose condition $(H_2)$ violates and $\theta _0/d_1\in (\mu _k,\mu _{k+1})$ for some $k\ge 1$, and $\sigma _k=\sum _{i=1}^{k} m(\mu _i)$ is odd, then there exists a positive constant $d^\star$, such that model (1.2) has at least one positive non-constant solution for all $d_2>d^\star$.

Proof

If $\theta _0>0$, then $(H_2)$ violates, which follows that for sufficiently large $d_2$ condition (4.9) holds and suggests $\mu _{+}(d_1,d_2)>\mu _{-}(d_1,d_2)$. Furthermore, we have

$$\begin{aligned} \lim _{d_2\rightarrow \infty } \mu _{+}(d_1,d_2)=\frac{\theta _0}{d_1},\ \ \lim _{d_2\rightarrow \infty } \mu _{-}(d_1,d_2)=0. \end{aligned}$$

For $\theta _0/d_1\in (\mu _k,\mu _{k+1})$, there exist $d_0>0$ which gives

$$\begin{aligned}&\mu _{+}(d_1,d_2)\in (\mu _k,\mu _{k+1}),\ \ \nonumber \\&0\le \mu _{-}(d_1,d_2)<\mu _1,\ \forall d_2>d_0. \end{aligned}$$

(4.10)

Theorem (8) states that for $d>d_0$ with $d_1=d$ and $d_2>d$ has no positive non-constant solution. If $d>0$ is large enough and $\frac{\theta _0}{d_1},\mu _1$ then there exists $d^\star >d$ such that

$$\begin{aligned} 0<\mu _{-}(d_1,d_2)<\mu _{+}(d_1,d_2)<\lambda _1 ,~~~~~\ \forall d_2>d^\star , \end{aligned}$$

(4.11)

Now we show that $d_2>d^\star$ ensures at least one positive non-constant steady state. To this end, we assume that the assertion is not true for some $d_2>d^\star$. We derive contradiction by homotopy argument. For $d_2=d_2^\star$ and $t\in [0,1]$, we define

$$\begin{aligned} D(t)=\begin{pmatrix} td_1+(1-t)d^\star &{} 0\\ 0 &{} dt_2+(1-t)d^\star \end{pmatrix} \end{aligned}$$

and considering

$$\begin{aligned} \left\{ \begin{aligned} -\varDelta \mathbf{u }&=D^{-1} \mathbf{F}(u) , ~~~~~~~~~~ x\in \varOmega ,\\ \frac{\partial \mathbf{u }}{\partial n}&=0 ~~~~~~~~~~~~~~~~~~~~~~ x\in \partial \varOmega . \end{aligned} \right. \end{aligned}$$

(4.12)

Here we note that u is a non-constant positive solution of (1.2) if and only if it is solution of (4.12) at $t=1$. For $t\in [0,1]$, u is a non-constant positive solution of (4.12) if and only if it is solution of the following

$$\begin{aligned} h(\mathbf{u },t):&=\mathbf{u }-(\mathbf{I }-\varDelta )^{-1}\{D^{-1}(t) \mathbf{F}(u) +\mathbf{u }\}\nonumber \\&=0, ~~~~ \text {on} \ ~~~ \mathbf{X }. \end{aligned}$$

(4.13)

It is important to note

$$\begin{aligned} h(\mathbf{u },1)=f(d_1,d_2;\mathbf{u }),~~~~h(\mathbf{u },0)=f(d,d^\star ;\mathbf{u }),\\ D_uf(d_1,d_2;E^\star )=\mathbf{I }-(\mathbf{I }-\varDelta )^{-1}(D^{-1}J+I)=0,\\ D_uf(d,d^\star ;E^\star )=\mathbf{I }-(\mathbf{I }-\varDelta )^{-1}({\tilde{D}}^{-1}J+I)=0, \end{aligned}$$

where ${\tilde{D}}=\begin{pmatrix} d &{}0\\ 0&{} d^\star \end{pmatrix}$ and f(., .; .) is defined in (4.8). From (4.10) and (4.11), we have

$$\begin{aligned} \varPi (d_1,d_2)\cap S_{p}&=\{\mu _0,\mu _1,\mu _2....\},~~ \text {and}\ \varPi (d,d^\star )\\&\quad \times \cap S_{p}=\emptyset . \end{aligned}$$

It is known that $\sigma _k$ is odd; thus, from Lemma (4),

$$\begin{aligned} \text {index}({h(.;1), E^\star })&=\text {index}({\bar{f}}(d_1,d_2;.),E^)\star \\&=(-1)^{\sigma _q}=-1,\\ \text {index}({h(.;0), E^\star })&=\text {index}({\bar{f}}(d,d^\star ;.),E^)\star \\&=(-1)^{0}=1. \end{aligned}$$

If $C^\star$ and ${\bar{C}}$ be positive constants satisfying $C^\star<u(x),v(x)<{\bar{C}}$ on $x\in {\bar{\varOmega }},\ t\in [0,1]$ where $C^\star =c^\star (d,d_1,d^\star ,d_2,\varLambda )$ and ${\bar{C}}=\bar{(}C)(d,d^\star , \varLambda )$. $\square$

We define

$$\begin{aligned} \varSigma =\{\mathbf{u }\in \ \mathbf{X }:C^\star<u(x),v(x)<{\bar{C}}, \ x\in {\bar{\varOmega }}\}. \end{aligned}$$

Then $h(\mathbf{u };t)\ne 0$ for all $u\in \partial \varSigma$. Now from Leray–Schauder degree, we have

$$\begin{aligned} \deg (h(.;0),\varSigma ,0)=\deg (h(.;1),\varSigma ,0) \end{aligned}$$

(4.14)

However, we note that $h(u;0)=0$ and $h(u;1)=1$ have a unique positive solution $E^\star$ in $\varSigma$. From the above, we obtain

$$\begin{aligned} deg(h(.;0),\varSigma ,0)&=\text {index}(h(.;0),E^\star )=1,\\ deg(h(.;1),\varSigma ,0)&=\text {index}(h(.;1),E^\star )=-1, \end{aligned}$$

which contradicts (4.14). This completes the proof.

4.4 Diffusion-driven instability

Now we investigate condition for occurrence of diffusion-driven instability, i.e., Turing instability for the system (1.2). Turing instability is a kind of diffusive instability which occurs due to the diffusion. So if homogeneous steady state of ODE system becomes unstable due to action of diffusion, then diffusion-driven instability occurs in the system. Recalling the hypothesis, conditions $\varPhi _1(k^2)>0$ or $\varPhi _1(k^2)<0$ insure diffusion-driven instability. Since we are interested in the growth of unstable patterns, we should select the largest perturbation $\lambda _k$ calculated in (4.5).

We observe that $\varPhi _2(k^2)$ defined in (4.4) is parabola and contains minimum value at $k_c^2$. It is important to determine the sign of $\varPhi _2(k^2)$ to examine the stability of the $E^\star$. Note that Turing instability may occur for some diffusion coefficients at which $\varPhi _2$ becomes negative for some critical positive wave number $k_c^2$. The following outcomes reveal that $\varPhi _2(k^2)>0$ is minimum at critical value $k_c^2$. Let $h=k^2$, now

$$\begin{aligned} \frac{\mathrm{d} \varPhi _2}{ \mathrm{d}k^2}&=2d_1d_2k^2+(d_1\theta _1-d_2\theta _0)=0\ \ \text {and}\ \ \frac{\mathrm{d}^2 \varPhi _2}{\mathrm{d}{k^2}^2}\\&=2d_2d_2>0. \end{aligned}$$

It is easy to calculate

$$\begin{aligned} k_c^2=\frac{d_2\theta _0-d_1\theta _1}{2d_2d_2} \end{aligned}$$

(4.15)

where $G(d_1,d_2)=d_2\theta _0-d_1\theta _1>0$. The minimum value of $\varPhi _2$ at $k_c^2$ is calculated as

$$\begin{aligned} \varPhi _2(k_c^2)&=-\frac{1}{4d_1d_2}\\&\quad \times \big [(d_1\theta _1-d_2\theta _0)^2+4d_1d_2(\theta _0\theta _1+\sigma s)\big ] \end{aligned}$$

Define $\eta =d_2/d_1$ and let

$$\begin{aligned} \digamma (d_1,d_2)&=(d_1\theta _1-d_2\theta _0)^2+4d_1d_2(\theta _0\theta _1+\sigma s)\\&=\theta _0^2 \eta ^2+2(\theta _0\theta _1+2s\sigma )\eta +\theta _1^2. \end{aligned}$$

Now

$$\begin{aligned} \digamma (d_1,d_2)=0&\Leftrightarrow \theta _0^2 \eta ^2+2\eta (\theta _0\theta _1+2s\sigma )+\theta _1^2=0, \end{aligned}$$

(4.16)

$$\begin{aligned} G(d_1,g_2)=0&\Leftrightarrow \eta =\frac{\theta _1}{\theta _0}\equiv \eta ^\star . \end{aligned}$$

(4.17)

Since $(\theta _0\theta _1+s\sigma )<0$ and $\sigma <0$. Thus, we have

$$\begin{aligned} 4(\theta _0\theta _1+2s\sigma )^2-4\theta _0^2 \theta _1^2=16s\sigma (\theta _0\theta _1+s\sigma )>0. \end{aligned}$$

Then $\digamma (d_1,d_2)=0$ has the following two positive roots

$$\begin{aligned} \mu _{1,2}=\frac{-(\theta _0\theta _1+ 2s\sigma )\pm 2s\sqrt{s\sigma (\theta _0\theta _1+s\sigma )}}{\theta _0^2} \end{aligned}$$

Therefore, when $d_2/d_1>\mu _1$, we have $\varPhi _2(k_c^2)<0$; thus, $E^\star$ is unstable. Hence, Turing instability occurs under the following condition:

$$\begin{aligned} \frac{d_2}{d_1}>\frac{-(\theta _0\theta _1+ 2s\sigma )+ 2s\sqrt{s\sigma (\theta _0\theta _1+s\sigma )}}{\theta _0^2}. \end{aligned}$$

(4.18)

then $\varPhi _2^2<0$. Thus, under condition (4.18) Turing instability occur at critical wave number $k_c^2$. If $k_\mathrm{max}$ refers maximum wave number, then wavelength ($\lambda$) can be calculated as $\lambda =\frac{2\pi }{k_\mathrm{max}}$ where $k_\mathrm{max}$ can be calculated from (4.15). The corresponding wave velocity can be obtained by $v=\frac{\delta x}{\delta _t}$ where $\delta x$ is the distance of the pattern passing in the time interval $\delta t$. The existence of diffusion-driven instability can be summarized by the following theorem

Theorem 10

If conditions ($H_1$), ($H_2$) and

$$\begin{aligned} \big (d_2\theta _0^2+d_1(\theta _0\theta _1+ 2s\sigma )\big )^2> 4\sigma s^3 d_1^2(\theta _0\theta _1+s\sigma ) \end{aligned}$$

hold. Then Turing instability occur near $E^\star$.

Remark 1

Now we show numerical illustration for existence of Turing instability near $E^\star$ for some threshold of diffusion coefficient. We numerically show that the $\varPhi _2(k^2)$ becomes negative when $k^2$ crosses some threshold $k_c^2$ and Turing instability appears. To this end, we set the following set of parameters

$$\begin{aligned} \left\{ \begin{aligned} \alpha _1&= 0.54,\ \alpha _2 = 0.12,\ \kappa = 0.1,\ \gamma _1 = 0.1, \ m = 0.27,\\ \gamma _2&= 0.1, \ c = 0.27;\beta = 0.1,\ d_1=d_2=1 \end{aligned} \right. \end{aligned}$$

(4.19)

that confers that solutions starting from a nonnegative initial data of the diffusive system (3.1) converge to constant steady state $E^\star =(0.841,\ 1.332)$. Figure 7a illustrates the stabilization of coexistence steady states $E^\star$ for the solutions starting from a positive initial data. It has been observe that setting $d_1=0.005$, $d_2=0.5$ and keeping all other parameters same as in (4.19) the numerical value of the expression $\varPhi _2(k^2)$ becomes negative for some positive wave number $k^2>0$. We refer readers to look at Fig. 6.

Fig. 6

Occurrence of Turing bifurcation for positive $k^2$

4.5 Hopf instability

In this subsection, we discuss the scenario when the system (1.2) is uniform in space and temporal symmetricity of the system breaks down and gives rise to oscillations that are uniform in space and periodic in time. In Theorem 5, we obtain Hopf bifurcation occurs at $b=b^\star$ at $E^\star$ of system (3.1). It can be seen that $\theta _0=H(u,v,b^\star )$; thus, for the simplicity we choose $\theta _0$ as bifurcation parameter. It is important to note that change parameter in b results symmetric change in bifurcation parameter $\theta _0$. It is known that characteristic polynomial 4.4 has purely imaginary roots if and only if $\varPhi _1(k^2)=0$ and $\varPhi _2(k^2)>0$. To this end, we have

$$\begin{aligned} \theta _0=\theta _0^\star :=\theta _1+(d_1+d_2)k^2. \end{aligned}$$

From the above observations, we have

$$\begin{aligned} \varPhi _2(k^2)\big |_{\theta _0^\star }=-\big [d_2k^2(k^2+2\theta _1)+(\theta _1^2+s\sigma )\big ]. \end{aligned}$$

From (4.16), we know that $\theta _1^2+s\sigma <0$. Hence, $\varPhi _2(k^2)\big |_{\theta _0^\star }>0$ for any $k\in \mathbb {N_+}$ if and only if $|d_2k^2(k^2+2\theta _1)|<|(\theta _1^2+s\sigma )|$ holds. It is easy to verify that if $\theta _0>\theta _0^\star$ for $k\in \mathbb {N_+}$, then positive steady state $E^\star$ is stable where as becomes unstable for $\theta <\theta _0^\star$. To verify transversality condition, let $i\omega$ is the purely imaginary root of (4.4). Then we have

$$\begin{aligned} -\omega ^2+\varPhi _2(k^2)=0,\ \varPhi _1(k^2)=(\theta _0-\theta _1)\\ \quad -(d_1+d_2)k^2=0. \end{aligned}$$

Then Hopf bifurcation occurs at $\theta _0=\theta ^\star$. Substituting $\lambda =p_1+iq_1$ in (4.4) and separating real and imaginary parts, we obtain

$$\begin{aligned} \begin{aligned} p_1^2-q_1^2-p_1\varPhi _1(k^2)+\phi _2(k^2)=0,\\ 2p_1q_1-q_1 \varPhi _1(k^2) =0. \end{aligned} \end{aligned}$$

(4.20)

Derivative of (4.20) along $\theta _0$ gives

$$\begin{aligned} \frac{\mathrm{d}p_1}{\mathrm{d}\theta _0}\big |_{\theta _0=\theta _0^\star }=-\frac{1}{2}<0. \end{aligned}$$

5 Taxis-driven instability

In this section, we extend the diffusive model (1.2) by take into account chemically mediated movement of predator toward the odor of the prey. We study the role of prey pheromones, chemical cues and chemical signaling on the dynamical behavior of the predator–prey model. Recently, Lee et al. [30] have studied effect of prey taxis on pattern formation accounting direct movement of predator toward the prey which stabilizes the system and known for annihilating Turing instability effect for large prey taxis sensitivity. Keeping this fact in our mind, we try to investigate how chemical signaling in predator–prey relation may alter the population dynamics. To this end, we assume that prey releases some chemical c(x, t) (e.g., pheromones) with rate $\gamma$ which vanishes with rate $\mu$. Taking into account the direct movement of predator toward the gradient of chemical released by prey, we try to examine the drivers that are responsible for patchiness in spatial distribution of species. We refer readers to look at papers [31,32,33] for the evidence of olfactory predators in the ecosystem. Recently, Tello and Wrzosek [34] have proposed a predator–prey model which accounts direct movement of predator toward the gradient of some chemical released by prey. This direct movement of predator toward the gradient of chemical released by prey is called indirect prey taxis. Ahn and Yoon [35] studied the global existence and uniform boundedness for a predator–prey model with general functional response and indirect prey taxis. Mishra and Wrzosek [36] proved global-in-time existence of solutions of an indirect taxis model with non-monotonic functional response and investigated the drivers responsible for pattern formation in predator–prey model. The revised predator model after incorporating indirect prey taxis can be written as

$$\begin{aligned} \left\{ \begin{aligned} \frac{\partial u }{\partial t}&=d_1 \varDelta u+ \frac{\alpha _1 u}{1+\kappa v}-\beta u^2-\frac{\gamma _1 uv}{bu^2+c}, ~~~~~~~~~~~x \in \varOmega ,\ t>0,\\ \frac{\partial v }{\partial t}&=d_2 \varDelta v-\chi \nabla \cdot (v\nabla c) +v\Big ( \alpha _2-\frac{\gamma _2 v}{m+u}\Big ), ~~~~x \in \varOmega ,\ t>0,\\ \frac{\partial c }{\partial t}&=d_3 \varDelta c-\mu c+\gamma u,~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ x \in \varOmega ,\ t>0, \end{aligned} \right. \end{aligned}$$

(5.1)

with nonnegative initial and Neumann boundary condition

$$\begin{aligned}&u(x,0)>0,\ v(x,0)>0,\ c(x,0)>0, ~~~~~~~~~~~~~~~ x\in \varOmega ,\ t=0,\\&\frac{\partial u}{\partial n}= \frac{\partial v}{\partial n}=\frac{\partial c}{\partial n} =0~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~x \in \partial \varOmega ,\ t>0.\\ \end{aligned}$$

We refer [35, 37] to authors interested in global-in-time existence and boundedness of solutions for such models. Here we examine linear stability of model (5.1). To this end, we linearize the system (5.1) near the coexistence steady state $e_t^\star =(u^\star ,v^\star ,\frac{\gamma u^\star }{\mu })$ where ($u^\star$,$v^\star$) are constant steady state of the model (1.2) defined in (3.2). The linearized system at the coexisting steady state $e_t^\star$ is written as

$$\begin{aligned} \frac{\partial }{\partial t}\begin{pmatrix} U\\ V\\ C \end{pmatrix} =\begin{pmatrix} d_1\varDelta +\theta _0 &{} \sigma &{} 0\\ s &{} d_2\varDelta +\theta _1 &{} -\chi v^\star \varDelta \\ 0 &{} \gamma &{} d_3\varDelta -\mu \end{pmatrix}\begin{pmatrix} U\\ V\\ C \end{pmatrix} \end{aligned}$$

where $\theta _0,\ \theta _1, s, \sigma$ are defined in (4.3). Let $\{h_n\}_{n=0}^{\infty }$ be the sequence of eigenvalues of operator $-\varDelta$ with homogeneous Neumann boundary conditions s.t $0=h_0<h_1\le h_2\le .....$ and

$$\begin{aligned} \lim _{i\rightarrow \infty } \ h_j=\infty . \end{aligned}$$

Then from Fourier expansion there exists $\phi _j$ such that

$$\begin{aligned} U&=\sum _{j=0}^{\infty } u_j \phi _j(x) e^{\lambda _jt},\ \ V=\sum _{j=0}^{\infty } v_j \phi _j(x) e^{\lambda _jt},\ \ \\ C&=\sum _{j=0}^{\infty } c_j \phi _j(x) e^{\lambda _jt}, \end{aligned}$$

where $u_j,\ v_j$ and $c_j$ are some constant. With the straightforward analysis, we obtain a Jacobian matrix as

$$\begin{aligned} J_c=\begin{pmatrix} \theta _0-d_1h_j &{} \sigma &{} 0\\ s &{} \theta _1-d_2h_j &{} \chi v^\star h_j\\ 0 &{} \gamma &{} -\mu -d_3h_j \end{pmatrix} \end{aligned}$$

(5.2)

Let $\lambda _j^1,\ \lambda _j^2$ and $\lambda _j^3$ be the eigenvalues of Jacobian matrix (5.2). Under the Routh–Hurtwitz criteria, the eigenvalues of Jacobian matrix (5.2). Thus, $E^\star _t$ is locally asymptotically stable if and only of for all $j\in {\mathbb {N}}$

$$\begin{aligned} \varphi _j^1&=\mu - \theta _0 - \theta _1 + (d_1+ d_2 + d_3)h_j,\\ \varphi _j^2&=-(\theta _0\theta _1 + \mu \theta _1 + s\sigma +\mu \theta _0) + (d_1 \mu +d_2\mu \\&\quad - d_1\theta _1- d_2\theta _0 - d_3\theta _0 - d_3\theta _1)h_j \\&\quad + (d_1d_2+ d_1d_3 + d_2d_3)h_j^2 + \chi \gamma v^\star h_j,\\&=\varphi _j^{2,1}+ \chi \gamma v^\star h_j,\\ \varphi _j^3&=\mu \theta _0\theta _1 - \mu s \sigma +h\Bigg ( d_1\mu \theta 1 + d_2 \mu \theta _0 + d_3 s\sigma \\&\quad - d_3\theta _0\theta _1\Bigg )+ h^2( d_1d_2\mu - d_1d_3\theta _1 - d_2d_3\theta _0)h_j^2\\&\quad + d_1d_2d_3h_j^3+ \chi \gamma v^\star h_j( \theta _0 - d_1 h_j),\\&= \varphi _j^{3,1}- \chi \gamma v^\star h_j( \theta _0 - d_1 h_j),\\&\varphi _j^1 \varphi _j^2-\varphi _j^3=\varphi _j^1\varphi _j^{2,1}-\varphi _j^{3,1}-\chi \gamma v^\star h_j M_j, \end{aligned}$$

where $\displaystyle M_j=\gamma v^\star h_j(\mu -\theta _1+d_2h_j+d_3h_j)>0$. Now we assume that $e^\star _t$ is locally asymptotically stable without taxis (i.e., $\chi =0$). Since $\chi$ appears linearly in $\varphi ^2_j$ and $\varphi ^3_j$, then from the stability conditions we have

$$\begin{aligned}&\varphi _j^1>0,\ \varphi _j^2=\varphi _j^{2,1}>0,\ \varphi _j^3=\varphi _j^{3,1}>0,\ \\&\varphi _j^1\varphi _j^{2}-\varphi _j^{3}=\varphi _j^1\varphi _j^{2,1}-\varphi _j^{3,1}>0. \end{aligned}$$

For $\chi >0$, it is easy to verify that conditions $\varphi _j^1>0,\ \varphi _j^2>0,\ \varphi _j^3>0$ holds. However, for the $\displaystyle \chi >\chi ^\star$ where

$$\begin{aligned} \chi ^\star = \min _{j\in {\mathbb {N}}_{+}}\Big \{ \frac{\varphi _j^1\varphi _j^{2,1}-\varphi _j^{3,1}}{M_j}\Big \} \end{aligned}$$

(5.3)

we have $\varphi _j^1 \varphi _j^2-\varphi _j^3<0$. Therefore, for a nonnegative $\chi ^\star$ the steady state $e_t^\star$ losses its stability at $\chi = \chi ^\star$ and becomes unstable for $\chi > \chi ^\star$. The above analysis can be summarized by the following theorem.

Theorem 11

The homogeneous steady state $e_t^\star$ of taxis-induced system (5.1) is locally stable if $\chi <\chi ^\star$. Steady state $e_t^\star$ looses its stability and becomes unstable at threshold of $\displaystyle \chi \ge \chi ^\star$.

Remark 2

It is noteworthy that a positive $\chi ^\star$ exists for which the stability criterion $\varphi _j^1 \varphi _j^{2} - \varphi _j^{3}>0$ may be violated and the Hopf-type instability might be possible if $\varphi _j^1 \varphi _j^{2} - \varphi _j^{3} <0$ holds for some $\chi ^\star > 0$. Please note that that Turing instability does not occur in Gauss-type indirect prey taxis systems for $\chi >0$ and inhomogeneous Hopf-type instability may exist. We refer authors to look at the paper [38] for detailed information on taxis-induced inhomogeneous Hopf bifurcation in the models of predator–prey interaction.

6 Numerical pattern formation

In this section, we present numerical patterning mechanism of a diffusive predator–prey system. We have obtained results for 1D case as well as 2D case and examined the impact of model parameters on pattern formations. With the help of numerical pattern formation, we discuss the biological and physical drivers responsible for the pattern formation. We use MATLAB PDEPE toolbox to obtain 1D numerical results, and two-dimensional model is simulated with the help of finite difference scheme. We take $\delta t=0.01$ and $\delta x=\delta y=0.1$ to increase the numerical accuracy of applied discretization scheme. For detail of the discretization scheme and finite difference method, we refer [39].

First we discuss existence of constant steady state and Hopf instability in 1D case. We take initial data of the form $u(x,y,0)=u^\star +0.01\cos (\frac{2\pi x}{50}), \ v(x,y.0)= v^\star +0.01\cos (\frac{2\pi x}{50})$ satisfying Neumann boundary condition and run simulations for the parameter set (4.19). We observe constant steady state for the set of parameter (4.19) and diffusion $d_1=d_2=1$. It is observed that constant steady becomes unstable at $\theta _0=\theta _0^\star$ via Hopf bifurcation and periodic solution is shown in Fig. 7b for $b=0.2$. It is important to note that these patterns are homogeneous in space and periodic in time. In the second panel of Fig. 7, we present space–time patterns for prey and predator, respectively. Figure 7a confirms locally stable homogeneous steady and its destabilization through Hopf bifurcation is presented in Fig. 7b.

Fig. 7

a Constant steady state at $b=0.1$, b periodic patterns for prey (u) and predator (v) induced by Hopf instability at $b=0.2$. Other parameters are same as in (4.19) except $d_1=d_2=1$

Fig. 8

Spatiotemporal pattern for prey (u), predator (v) and chemical (c) with $\chi =0.97>\chi ^\star$ and $\mu =0.2,\ \gamma =1,\ d_1=1, d_2=d_3=0.001$ for the taxis system (5.1); and other parameters same as in (4.19)

Now we present interesting numerical results for the taxis-induced system (5.1). We take symmetric initial data of the form $u_0=0.4685+0.01\cos (\frac{4\pi x}{50}),\ v_0= 0.95+0.01\cos (\frac{4\pi x}{50}),\ c_0= 2.45+0.01\cos (\frac{4\pi x}{50})$that confirms the stability of homogeneous steady state at $\chi =0$ (i.e., without taxis). We computed the critical chemical sensitivity $\chi ^\star$ for the a particular set of parameter value. From Fig. 8, it is observed that the steady state solution is locally stable for the $\chi <\chi ^\star$ and looses its stability and spatiotemporal pattern occurs for $\chi >\chi ^\star$. It is important to notice that spatiotemporal pattern is inhomogeneous in space and periodic in time. One can see formation of various patches of prey, predator and chemical in space and it is interested to see that similar patches evolve after some period of time.

Fig. 9

Evolution of target wave patterns for prey (u) and predator (v) in time a $t=0$ b $t=500$, c $t=2000$, d $t=6000$. Other parameters are $\alpha _1 = 0.54,\ \alpha _2 = 0.12, \kappa =0.1, \ b = 0.1,\ \gamma _1 = 0.1,\ \gamma _2 = 0.1,\ m = 0.27, \ c = 0.27;\beta = 0.1, d_1=0.005,\ d_2=0.5$

6.1 Turing pattern formation

Now we present some interested Turing patterns observed in 2D square domain $\varOmega \in (50,\ 50)$. In our previous discussion, we have shown occurrence of Turing instability (see Fig. 6). In this section, we run simulations for a large time and take parameters in a way that guarantees the Turing pattern. Our main objective of numerical simulations is to show the transition of pattern based on biological as well as physical factors. As we are interested in the drivers of pattern formation, we have thus run simulations to examine the effects of ecological parameters and the prevalence of individuals in Turing domain.

6.1.1 Emergence of target wave pattern

We present numerical simulations obtained from initial data

$$\begin{aligned} \left\{ \begin{aligned} u(0)&=u^\star +0.1*e^{-(x-30)^2-(y-30)^2},\\ v(0)&=v^\star +0.1*e^{-(x-30)^2-(y-30)^2} \end{aligned} \right. \end{aligned}$$

(6.1)

to examine the evolution of Turing patterns in the advancement of time. It is worth mentioning that numerical results for target wave patterning are obtained for the parameters set which insures pure Turing patterns. To this end, we run simulations for various time intervals in the Turing domain and presents snapshots for prey and predator distribution over square domain $\varOmega$. The first snapshot presented in Fig. 9a is captured at time $t=0$ which shows the shape of initial data (6.1) for the simulations. In Fig. 9b, we observe formation of very beautiful concentric ring-shaped structure at the center of the domain for time $t=500$. It is worth mentioning that red color represents the high density of individuals and the blue shows low concentration of individuals. The color of the snapshot suggests that both prey and predator species form ring-shaped patches of high density. In the advancement of time $t=2000$, we observe formation of various concentric rings across the domain (see fig. 9c). A keen observation on this figure suggests that prey is able to form packed ring-shaped flocks of very high density and remains very few outside the ring. However, predator’s density at prey dominant area remains low. As the time progresses and at reaches $t=6000$, individuals of prey and predator prefer to stay in ring-shaped packed group (see Fig. 9d). Obviously the width of ring is getting wider and emergence of dark blue symmetric structure outside the boundary of ring emerges for both prey and predator’s distribution. It is worth noticing that all emerged concentric target waves in Fig. 9d have same center. It is worth mentioning that we observe no change in patterns for a larger time after $t=6000$. Therefore, we can comment that a solution starting from initial data (6.1) showing a patch both prey and predator at the middle of a square domain forms stationary concentric target wave-shaped distribution. It is very interesting to see that prey concentration on target wave boundary remains very high. However, predator also try to settle down with prey on the boundary but declines over time. The formation of target wave patterns suggests that prey individuals move toward the nearest target wave and forms more packed group. In the meantime, predator also follows the prey and tries to aggregate the herd of prey for suggestive predation. The target waves in ecology are driven by pursuit and invasion and demonstrate the migration between subpopulations which is a potential ecological significance.

6.1.2 Emergence of spot and stripe patterns

In order to obtain spot- and stripe-like patterns in Turing domain, we start numerical simulations from random initial data

$$\begin{aligned} u(0)=u^\star +0.1\xi ^{\star }\ v(0)=v^\star +0.1\xi ^{\star }, \end{aligned}$$

(6.2)

where $\xi ^{\star }\in (0,1)$ is any random number. Here we present snapshots captured at different time intervals and tried to show that solutions starting from random initial data (6.2) give rise to time-independent strip- and spot-like structure. In Fig. 10a, we present the structure of initial data (6.2) that represent the distribution of prey and predator individuals at time $t=0$ over the square domain $\varOmega \in [0,50]\times [0,50]$. The snapshot captured at time $t=50$ suggests that randomly distributed prey individual quickly gathers and forms their packed groups. Red spots in Fig. 10b show the flocks of individuals, and blue color represents low density of individual. For a larger time $t=500$, one can see emergence of beautiful pattern that contain mixture of stripes and spots. However, it is noteworthy that stripe-like structure dominates in the square domain (see Fig. 10c). We observe that stripes are packed with very high density of prey individuals and predator also settles down inside the same stripe but predator’s average density remains comparatively low inside the stripe structure. Figure 10d is captured at the time $t=3000$ and it turns out that these pattern are steadily with time and temporally independent. We observe that patterns shown in Fig. 10d are steady state patterns obtained in Turing domain. It is worth notifying that hot stripes dominate in the domain. Stripe and spot patterning is very ubiquitous in nature and play role in fish patterning.

Fig. 10

Evolution of Turing patterns for prey (u) and predator (v) in time a $t=0$ b $t=50$, c $t=500$, d $t=3000$ e $t=6000$. Other parameters are $\alpha _1 = 0.54,\ \alpha _2 = 0.12, \kappa =0.1, \ b = 0.1,\ \gamma _1 = 0.1,\ \gamma _2 = 0.1,\ m = 0.27, \ c = 0.27;\beta = 0.1, d_1=0.005,\ d_2=0.5$

Now we try to investigate the impact of diffusion coefficients on the pattern formations. To this end, we start numerical simulations with the randomize initial data (6.2) and keeping all other parameters same as in Fig. 10 except $d_1=0.007$ and changes in the value of $d_2$. We run simulations for a large time period $t = 10000$ which guarantees the non-dependence of the patterns on time giving us a pure diffusion-driven Turing patterns. In Fig. 11a, we observe the emergence of spot patterns of both prey and predator individuals for predator’s diffusion coefficient $d_2=2$. It is very interesting to see formation of red spots for both the predator and prey individuals. This result biologically reveals that prey and predator individuals settle down at the same location where prey density dominates over predator. Increasing the diffusivity rate of the predator suggests that spot-like structure transits into spot–stripe mixture patterns. One can see mixtures of stripe- and spot-like structures in Fig. 11b obtained for the diffusive rate $d_2=0.35$. In our case, both predator and prey individuals settle down at the same location and forms a stripe pattern. However, few hot spots can be also seen. In Fig. 11c, we observe stripe pattern obtained for $d_2=0.55$. It has been observed that as we increase the diffusion rate of predator, mixtures of stripes and spots transits into pure stripe patterns. It is worth mentioning that red stripe represents a location at which individuals settle down and movement of individuals inside stripe structure is possible. Therefore, the transition from the pattern of spots to the striped pattern is expected as the predator prevalence rate increases.

Fig. 11

Snapshots of Turing patterns for prey (u) and predator (v) at time $t=10000$ for parameter set $\alpha _1 = 0.54,\ \alpha _2 = 0.12, \kappa =0.1, \ b = 0.1,\ \gamma _1 = 0.1,\ \gamma _2 = 0.1,\ m = 0.27, \ c = 0.27;\beta = 0.1, d_1=0.007$ except a $d_2=0.2$ b $d_2=0.35$, c $d_2=0.55$

Fig. 12

Snapshots of Turing patterns for prey (u) and predator (v) at a $\kappa =0.05$ b $\kappa =0.1$, c $\kappa =0.15$. Other parameters are $\alpha _1 = 0.54,\ \alpha _2 = 0.12,\ b = 0.1,\ \gamma _1 = 0.1,\ \gamma _2 = 0.1,\ m = 0.27, \ c = 0.27;\beta = 0.1, d_1=0.007,\ d_2=0.3$

Fig. 13

Snapshots for Turing patterns of prey (u) and predator (v) at a $b=0.05$ b $b=0.1$, c $b=0.15$. Other parameters are $\alpha _1 = 0.54,\ \alpha _2 = 0.12,\ \kappa = 0.1,\ \gamma _1 = 0.1,\ \gamma _2 = 0.1,\ m = 0.27, \ c = 0.1;\beta = 0.1, d_1=0.007,\ d_2=0.3$

6.1.3 Role of fear and grouping on patterns

Now we present snapshots to show how fear and inducible defenses alter the population dynamics in a diffusive predator–prey system. We present snapshots of time-independent solutions that are emerging due to changes in the respective parameters. We take random initial data of the form (6.2) and run simulations for large time $t=10000$. First we present snapshots observed for various values of $\kappa =0.05, 0.1 and \ 0.15$ in Fig. 12 to show impact of fear on the density distribution of populations. In Fig. 12a, we observe emergence of hot labyrinth prey patterns for $\kappa =0.05$ in the square domain $\varOmega \in [0,50] \times [0,50]$. It is important to note that labyrinth evolves in time continuation and are spreading labyrinth. We observe that prey is highly distributed inside the structure, which suggests that these hot streaks develop due to aggregation. It has been observed that predators try to follow prey and forms very similar complex structure. In the next simulation, we take $\kappa =0.1$ and keep other parameters same. We observe that both prey and predator pattern transit into stripe–spot mix pattern (see Fig. 12b). This result informs that prey are more disciplined when intensity of fear is high enough. We observe formation of stripes into hot spots for $\kappa =0.15$. However, a few stripes are shown in Fig. 12c. This transition of patterns depending upon the parameter $\kappa$ suggests that fear has immense impact on the distribution of species which alters the population dynamics of both prey and predator individuals. Now we present simulations for different values of parameter b in Turing domain and try to capture snapshots at various time of intervals. All the parameters and initial data are kept same as in Fig. 12. Figure 13a is obtained for the parameter $b=0.05$. We observe very small population density of prey, whereas predator population dominates in the space. For a higher value of parameter $b=0.1$, we observe formation of hot spot and stripe patterns. It is important to note that prey forms these structures and the color bar suggests high population density inside these structures (see Fig. 13b). A very similar structure is obtained for the predator’s population. In Fig. 13b, we observe formation of very complex structure when the parameter b is high enough. It is worth mentioning that we observe emergence of complex labyrinth structure for the parameter $b=0.12$ which is very close to the bifurcation threshold. We observe the transition of mixtures of stripe–spot patterns to complex labyrinth that means as defense efficiency b increases the prey individuals form group. A combination of hot as well as cold spot can be seen in distribution of predator population. In case of more dangerous prey, we observe that prey form a symmetrical cold stripe-like flock pattern which is surrounded by red color that represents high density of prey. This snapshot biologically reveals that somehow prey individuals prefer to stay in flocks that helps them to fight against predator and gives a chance for the survival. It is interesting to note that predator dominates in domain except the area covered by highly distributed prey (see Fig. 13c).

6.2 Taxis-induced patterns

In this subsection, we present 2D patterns induced by taxis-driven instability for extended chemo-attraction model (5.1) in Fig. 14. We use FreeFem++ package in order to obtain accuracy in numerical results we run simulations for the small time and space steps $\delta t=0.01,\ \delta x=\delta y=0.1$. Figure 14 shows that taxis-driven instability gives rise to spiky solutions which biologically represent patch of population in the domain. It has been observed that if the chemo-attractive sensitivity is high enough and both prey and predator dominate initially at the middle of the domain then interior spike is evolved. Figure 14a is obtained for initial data $u(0)=u^\star +e^{-(x-5)^2-(y-5)^2},\ v(0)=u^\star , c(0)=v^\star$ and parameter set are taken as in (4.19) except $\chi = 2.0,\ \mu = 0.2,\ \gamma = 1, \ d_1 = 1,\ d_2 = d_3 = 0.001$. Please note that taking different initial data has biggest impact on the taxis-induced patterns. It has been observed that boundary spikes emerge for initial date of the form $u(0)= u^\star + 0.1 \cos (\frac{\pi x}{10}) \cos (\frac{\pi x}{10}), v(0) = v^\star , \ c(0) = c^\star$. These figures suggest that spiky periodic solutions are possible when predator–prey model is influence by indirect prey taxis. However, it is worth mentioning that these patterns are space–time pattern.

Fig. 14

Spatiotemporal spiky patterns for prey (u) and predator v a Emergence of taxis-driven interior spikes, b emergence of taxis-driven boundary spikes in two-dimensional chemo-attraction-induced model (5.1)

7 Conclusion

In natural populations, fear of predation is evident and in response the prey exhibits some form of self-defense to protect themselves. Both predator and prey species have adaptations—beneficial characteristics generated by natural selection—that help them perform better in their roles. For example fear response is widely detected in small fish and schooling is one of the primary ways many fishes adopt to defend themselves from predators [14]. Motivated by this fact, we have proposed and studied impact of prey schooling on fearful prey in a diffusive predator–prey model. We have incorporated prey growth function subjected to improvements through introduction of fear effects. We particularly focused on realistic biological phenomena that may alter the population dynamics. We have incorporated simplified Monod–Haldane functional response that resembles non-monotonic impact of prey schooling on predation and studied dynamical properties of model. The dynamics of proposed Leslie–Gower-type model has been investigated at various population states. We have derived the stability conditions and established restrictions for occurrence of Hopf bifurcation for temporal model. The direction of Hopf bifurcation is verified by calculation of Lyapunov coefficient at the bifurcation parameter.

In case of diffusive system, we have shown that system is well posed and local stability of steady states is investigated with the help of linearization technique. Linear stability conditions confer that steady state remains locally stable if predator diffuses slowly and local temporal system is locally stable. We obtain conditions that ensure the existence of non-constant steady state at some threshold of parameters. It is also observed that coexistence steady state of diffusion-induced system destabilizes and Hopf instability appears at some threshold of defense intensity. Turing-like instability is obtained under certain condition but with the restriction that predator should be very quicker than prey individuals. It is worth notifying that not only Turing- and Hopf-like instabilities give rise to spatiotemporal pattern, but also indirect prey taxis gives rise to spatiotemporal patterns. Our analysis demonstrates that chemical clues play very important role on pattern formations when predator makes taxis toward the chemical released by the prey. However, direct taxis of predator toward the gradient of prey promotes stability in predator–prey relations and pattern formation is unobtainable for high taxis sensitivity [30]. On the other hand, indirect prey taxis in Gauss-type predator–prey models also annihilates Turing spatial patterns and gives rise to spatially inhomogeneous Hopf bifurcations [38]. Numerically we observe stable attractor and a stable limit cycle for the non-diffusive version of the system. Sensitivity analysis confirms that fear is able to stabilize the coexisting system; however, group defense destabilizes the predator–prey model. There are few studies that supports our results [16, 20]. However, Pat et al. [40] investigated destabilizing effect of fear in a predator–prey model with mutual interference among predators. Destabilization effect of prey fear is also studied with strong adaptation of adult prey and a special class of predator–prey models that constitute prey refuge [18, 41]. Various Turing patterns such as spot, stripes and labyrinth are observed without indirect taxis. Our most important finding is emergence of target wave patterns that reveals the applicability of this study in the real-world system. Similar target wave pattern is seen in special class of reaction diffusion prey–predator models [42]. Periodic travelling wave mechanism is observed in real-world field studies of vole and Canadian lynx populations; and in special class of cyclic predator–prey models [43, 44]. We numerically observed inhomogeneous space–time patterns when indirect prey taxis is presented in predator–prey system. However, it remains unknown whether non-constant time-independent pattern are possible or not in this case which may me a motivation for the our readers in this direction for further investigation. Our numerical simulations demonstrate that physical drivers such as initial concentrations, diffusion rates and chemical mediated signaling play important role in the shaping of spatial patterns and are also able to alter the spatial distribution of species in any ecosystem.

Finally this study reveals that biological drivers such as anti-predation, fear response and invasion scenario play very important roles in predator–prey interaction and responsible for altering the population dynamics. However, physical drivers such as movement of individuals are responsible for pattern formation. This study also demonstrates that fear-induced anti-predation helps in prey survival, although fear predisposition is characteristic to restrict the development of prey individuals. One of the main result of this paper states that the fear response and indirect prey taxis is capable of contrasting the diffusion-driven instability effect of group defense in predator–prey systems. Biological implication of this study reveals that staying in schools is always beneficial and the aspect of the “safety in numbers” theory revolves around the fact that even if a predator attacks a school, the odds are low that any one prey will be the one the predator captures but predator has to bear its consequences.