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

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.


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 antipredation 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 (Hel-ogale 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 feardriven 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  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 τ = t h + t f u. Thus, Holling disk equation can be rewritten as f (u) = au 1 + at h u + at f u 2 , (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) f (u) = 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 nonmonotonic 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].
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 preda-tor 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 schoolingdependent 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 predatorprey 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 predatorprey 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 In model (1.2), the term g(u, v, κ) = α 1 u 1 + κv represents growth rate of prey as a function of preda-tor density. Here g(g, v, 0) = α 1 u, g(u, 0, κ) = α 1 u, lim κ→∞ g(u, v, κ) = 0, lim v→∞ g(u, v, κ)=0, ∂g ∂κ < 0 and ∂g ∂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 α 1 measures the growth rate of prey and nonnegative parameter κ represents the intensity of fear. The prey growth function g(u, v, κ) without fear (i.e., κ = 0) is g(u, v) = α 1 u, which represents exponential growth in prey. We incorporate non-monotonic prey capture response to model group defense in predatorprey system with encounter rate γ 1 . The parameter β measures mortality rate due to intraspecific competition, b measures intensity of defense, c is saturation constant and α 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 − γ 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 lim u→∞ α 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-intime 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.

Existence of global positive solutions
Let us denote x + : x∨0, x − : (−x)∨0 for x, y ∈ IR and V + = {v ∈ V : v ≥ 0} which is positive cone of a partial order vector space (V, ≥). Assume Z = (u, v) ∈ X and its norm is defined as The system (1.2) can be written as the following abstract equation: holds Duhamel formula for t ∈ [0, t max )

4)
where operator A p represents a closure of A ∈ X. Moreover, if t max < ∞ then 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 t A p . Furthermore, for t > 0, we have Moreover, here we observe that B : Λ → X is locally Lipschitz on a bounded set. Thus, from Theorem 7.2.1 in [26], we can easily obtain desired result.

Lemma 2
For nonnegative initial data and for p > N , , v(x, t)) of the system (1.2) are also nonnegative.
Proof The system (1.2) in positive cone can be written as with initial and boundary conditions We multiplyū − andv − , respectively, in (2.5) and applying by part integration on domain Ω we obtain Therefore, for t ∈ (0, t max ), we obtain is the solution of the system (1.2) and Lemma (2) and we obtain u( . 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
Then lim t→inf p(x, t) = α for any x ∈ Ω.

Theorem 1
The solution (u(x, t), v(x, t)) of the system (1.2) satisfies the following inequalities Proof From the equation of prey growth in (1.2), From comparison inequality and Lemma (3), we can easily obtain first inequality. Let ∈ (0, 1] is sufficiently small arbitrary constant then there exists a T 1 > 0 such that u(x, t) ≤ α 1 β + for any x ∈ Ω and t > T 1 . Now for the second equation, we have From the comparison inequality and Lemma (3), for any ∈ (0, 1], there exists at From the arbitraryness of , we get lim sup This completes the proof. 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 (u(x, t), v(x, t) such that

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 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.

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 = ( α 1 β , 0) exists for positive values of α 1 and β.
Now we discuss the existence of coexisting equilibrium state E = (u , v ). From the second equation Substituting v in first equation of ODE system, we get polynomial of the form where 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 . Hence, unique coexisting equilibrium E = (u , v ) exists under some certain restrictions. We numerically choose a set of parameters which gives a unique interior equilibrium E . 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 Here both eigenvalues are positive. Hence, trivial equilibrium E 0 is always unstable. (ii) The Jacobian matrix calculated at E 1 = ( α 1 β , 0) Jacobian matrix J 1 has two eigenvalues −α 1 and α 2 of the opposite signs. Hence, predator-free equilibrium E 1 is a saddle point.
Here −α 2 and the Jacobian matrix J 2 reduces to 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 T 2 at α 1 = α 1 , respectively, then we obtain V = Again in order to confirm transcritical bifurcation we calculate Similarly following [27], one can (iv) Jacobian matrix at coexisting equilibrium point E is calculated as , The characteristic polynomial of J E is eigenvalues for the Jacobian matrix J E . Then E is locally asymptotically stable.
From the above analysis, we have the following theorem: Theorem 3 The ODE system has the following stability properties then coexisting equilibrium E is locally asymptotically stable.
Now we present the conditions for global stability E with the help of Lyapunov direct method. Global stability of solutions can be summarized by the following theorem:

Theorem 4 Suppose E is locally asymptotically stable and condition (A1) holds, then coexisting equilibrium E of ODE system (3.1) is globally asymptotically stable.
Proof We apply Lyapunov direct method to examine the global stability of coexisting equilibrium E of ODE system (3.1). Let us consider the Lyapunov function of the form Derivative along time t gives Here σ is a positive parameter to be chosen suitably later. Some mathematical calculations give Writing the above equation in the form of quadratics, we get . Therefore, coexisting equilibrium E of ODE system (3.1) is globally asymptotically stable initiating in positive invariant set Λ. This proves our result.

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 . Now we examine the Hopf bifurcation near the coexisting equilibrium E . In addition, we also discuss the direction of Hopf bifurcation at coexisting E when bifurcation parameter crosses its Theorem 5 E undergoes for a Hopf bifurcation when parameter crosses its critical value. b = b . Moreover, Hopf bifurcation is supercritical or subcritical depending upon ρ < 0 or ρ > 0, respectively.
Proof The Jacobian matrix at E is presented in Eq.
(3.6). When Jacobian matrix J E has pure imaginary eigenvalue of the form λ = ±ω and transversality condition holds, then Hopf bifurcation occur for b = b . Obviously, it is sufficient to show that Hopf bifurcation occurs for some critical value Setting φ 1 = 0 and performing some straightforward calculations, give the following threshold value for which conditions (i) and (ii) hold. Now we focus on the direction of Hopf bifurcation and stability of periodic solutions staring nearby E . In order to examine the stability of periodic solution, we calculate the Lyapunov coefficient at α 1 = α 1 . We translate E 0 to the origin by transformation u = u − p and v = v − q. Substituting these transformation in ODE model (3.1) and expanding with the help of Taylor series, we get and a i j , b i j are defined as The first Lyapunov coefficient defined in [28] is given by For system (3.1), we have 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 The stability of periodic solution depends upon Lyapunov coefficient ρ. The limit cycle is stable if ρ < 0 and unstable if ρ > 0. Therefore, Hopf bifurcation is supercritical if ρ < 0 and supercritical if ρ > 0.

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  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 = 0.165. In Fig. 4, we observe that solutions of u(t) and v(t) oscillated around E . This represents occurrence of limit cycle and destabilization of coexisting equilibria. Thus, Hopf bifurcation is verified for b > b . 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 and demonstrates destabilization of coexistence population  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. Now in order to determine the stability and direction of Hopf bifurcation, we calculate coefficients a i j and b i j to obtain the sign of Lyapunov coefficient ρ. Numerical simulation gives the following at b = 0.164 We obtain ρ = −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.

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 diffusiondriven 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).

Stability of steady state
In this section, we derive the Turing bifurcation conditions around the homogeneous equilibrium E 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 where k is wave number, λ k represents change in perturbation in time t and c k 2 and c k 2 are the amplitudes. Substituting the above transformations in Eq. (4.2), we get the following linearized form (4.1) The standard linear form of the PDE system (1.2) can be written as Here D = diag(d 1 , d 2 ) and L represents a linear oper- J represents Jacobian matrix at steady state E and defined as where θ 0 = A 11 , σ = A 12 , s = A 21 and θ = −A 22 are defined in (3.5). From the above Jacobian matrix, we get The characteristic polynomial of the above Jacobian matrix is where The roots of (4.4) are If (H 3 ) holds then θ 1 < 0. One can easily verify that σ < 0, s > 0 and θ 1 > 0. Then Φ 1 (k 2 ) < 0 and Φ 2 (k 2 ) > 0, we conclude that the characteristic polynomial (4.4) have roots with negative real part for all k ≥ 0. Hence, if d 1 d 2 > θ 0 θ 1 i.e predator is much faster than prey and conditions (H 3 ) and (H 4 ) holds then Φ 2 (k 2 ) is positive. Therefore, the coexistence steady state E with diffusion is locally stable. The above analysis can be summarized by the following theorem: 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:

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 Proof Let us assume that (u, v) is a positive solution of (1.2). Letφ = 1 |Ω| Ω ϕdx for ϕ ∈ L 1 (Ω). We multiply u −ū to the first equation of (1.2) and then integrating on Ω Here c 1 , are positive constants and ξ ∈ (u,ū) and μ ∈ (v,v). Similarly we have Adding the above quantities, we get for positive constant nu and some arbitrary constant > 0. Applying Poincáre inequality, we obtain Since μ 1 d 2 > α 2 , we may choose sufficiently small 4.3 Existence of positive non-constant steady states In this subsection, we discuss the existence of nonconstant steady state solutions depending upon diffusivity of individuals. For the simplicity, we denote 3), and f and g are defined in (4.6). The system (1.2) can be written as where where (I − Δ) −1 is inverse of (I − Δ) with homogeneous no-flux boundary condition. Straightforward calculations give For ∀X i if and only of λ is eigenvalue of the following matrix thenḡ(d 1 , d 2 , λ) = 0 has the following real eigenvalues Let us define where S p is positive spectrum of Δ on Ω with homogeneous NBC. We apply the following lemma to calculate In particular, ifḡ(d 1 , d 2 ; μ i ) > 0 for all μ i ≥ 0, then σ = 0. Proof If θ 0 > 0, then (H 2 ) violates, which follows that for sufficiently large d 2 condition (4.9) holds and suggests μ + (d 1 , d 2 ) > μ − (d 1 , d 2 ). Furthermore, we have For θ 0 /d 1 ∈ (μ k , μ k+1 ), there exist d 0 > 0 which gives (4.10) Theorem (8) (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. It is important to note ; u), It is known that σ k is odd; thus, from Lemma (4), We define which contradicts (4.14). This completes the proof.

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 diffusiondriven instability occurs in the system. Recalling the hypothesis, conditions Φ 1 (k 2 ) > 0 or Φ 1 (k 2 ) < 0 insure diffusion-driven instability. Since we are interested in the growth of unstable patterns, we should select the largest perturbation λ k calculated in (4.5). We observe that Φ 2 (k 2 ) defined in (4.4) is parabola and contains minimum value at k 2 c . It is important to determine the sign of Φ 2 (k 2 ) to examine the stability of the E . Note that Turing instability may occur for some diffusion coefficients at which Φ 2 becomes negative for some critical positive wave number k 2 c . The following outcomes reveal that Φ 2 (k 2 ) > 0 is minimum at critical value k 2 c . Let h = k 2 , now It is easy to calculate Since (θ 0 θ 1 + sσ ) < 0 and σ < 0. Thus, we have Then (d 1 , d 2 ) = 0 has the following two positive roots Therefore, when d 2 /d 1 > μ 1 , we have Φ 2 (k 2 c ) < 0; thus, E is unstable. Hence, Turing instability occurs under the following condition: then Φ 2 2 < 0. Thus, under condition (4.18) Turing instability occur at critical wave number k 2 c . If k max refers maximum wave number, then wavelength (λ) can be calculated as λ = 2π k max where k max can be calculated from (4.15). The corresponding wave velocity can be obtained by v = δx δ t where δx is the distance of the pattern passing in the time interval δt. The existence of diffusion-driven instability can be summarized by the following theorem Theorem 10 If conditions (H 1 ), (H 2 ) and hold. Then Turing instability occur near E .
Remark 1 Now we show numerical illustration for existence of Turing instability near E for some threshold of diffusion coefficient. We numerically show that the Φ 2 (k 2 ) becomes negative when k 2 crosses some threshold k 2 c and Turing instability appears. To this end, we set the following set of parameters that confers that solutions starting from a nonnegative initial data of the diffusive system (3.1) converge to constant steady state E = (0.841, 1.332). Figure 7a illustrates the stabilization of coexistence steady states E 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 Φ 2 (k 2 ) becomes negative for some positive wave number k 2 > 0. We refer readers to look at Fig. 6.

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 at E of system (3.1). It can be seen that θ 0 = H (u, v, b ); thus, for the simplicity we choose θ 0 as bifurcation parameter. It is important to note that From the above observations, we have From (4.16), we know that θ 2 1 + sσ < 0. Hence, Φ 2 (k 2 ) θ 0 > 0 for any k ∈ N + if and only if |d 2 k 2 (k 2 +2θ 1 )| < |(θ 2 1 +sσ )| holds. It is easy to verify that if θ 0 > θ 0 for k ∈ N + , then positive steady state E is stable where as becomes unstable for θ < θ 0 . To verify transversality condition, let iω is the purely imaginary root of (4.4). Then we have Then Hopf bifurcation occurs at θ 0 = θ . Substituting λ = p 1 +iq 1 in (4.4) and separating real and imaginary parts, we obtain Derivative of (4.20) along θ 0 gives

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 γ which vanishes with rate μ. 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 with nonnegative initial and Neumann boundary condition We refer [35,37] to authors interested in global-intime 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 = (u , v , γ u μ ) where (u ,v ) are constant steady state of the model (1.2) defined in (3.2). The linearized system at the coexisting steady state e t is written as Then from Fourier expansion there exists φ j such that where u j , v j and c j are some constant. With the straightforward analysis, we obtain a Jacobian matrix as Let λ 1 j , λ 2 j and λ 3 j be the eigenvalues of Jacobian matrix (5.2). Under the Routh-Hurtwitz criteria, the eigenvalues of Jacobian matrix (5.2). Thus, E t is locally asymptotically stable if and only of for all j ∈ N Now we assume that e t is locally asymptotically stable without taxis (i.e., χ = 0). Since χ appears linearly in ϕ 2 j and ϕ 3 j , then from the stability conditions we have For χ > 0, it is easy to verify that conditions ϕ 1 j > 0, ϕ 2 j > 0, ϕ 3 j > 0 holds. However, for the χ > χ where we have ϕ 1 j ϕ 2 j − ϕ 3 j < 0. Therefore, for a nonnegative χ the steady state e t losses its stability at χ = χ and becomes unstable for χ > χ . The above analysis can be summarized by the following theorem.

Theorem 11 The homogeneous steady state e t of taxis-induced system (5.1) is locally stable if χ < χ . Steady state e t looses its stability and becomes unstable at threshold of χ ≥ χ .
Remark 2 It is noteworthy that a positive χ exists for which the stability criterion ϕ 1 j ϕ 2 j −ϕ 3 j > 0 may be violated and the Hopf-type instability might be possible if ϕ 1 j ϕ 2 j − ϕ 3 j < 0 holds for some χ > 0. Please note that that Turing instability does not occur in Gauss-type indirect prey taxis systems for χ > 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.

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 MAT-LAB PDEPE toolbox to obtain 1D numerical results, and two-dimensional model is simulated with the help of finite difference scheme. We take δt = 0.01 and δx = δ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 +0.01 cos( 2π x 50 ), v(x, y.0) = v + 0.01 cos( 2π 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 θ 0 = θ 0 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. 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( 4π x 50 ), v 0 = 0.95 + 0.01 cos( 4π x 50 ), c 0 = 2.45 + 0.01 cos( 4π x 50 )that confirms the stability of homogeneous steady state at χ = 0 (i.e., without taxis). We computed the critical chemical sensitivity χ for the a particular set of parameter value. From Fig. 8, it is observed that the steady state solution is locally stable for the χ < χ and looses its stability and spatiotemporal pattern occurs for χ > χ . 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.

Turing pattern formation
Now we present some interested Turing patterns observed in 2D square domain Ω ∈ (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,

Emergence of target wave pattern
We present numerical simulations obtained from initial data 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 Ω. 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 ringshaped 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.

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 where ξ ∈ (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 Ω ∈ [0, 50] × [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. 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.

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 κ = 0.05, 0.1and 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 κ = 0.05 in the square domain Ω ∈ [0, 50] × [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 sug- gests 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 κ = 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 κ = 0.15. However, a few stripes are shown in Fig. 12c. This transition of patterns depending upon the parameter κ 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).

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 δt = 0.01, δx = δ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 chemoattractive 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 + e −(x−5) 2 −(y−5) 2 , v(0) = u , c(0) = v and parameter set are taken as in (4.19) except χ = 2.0, μ = 0.2, γ = 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 + 0.1 cos( π x 10 ) cos( π x 10 ), v(0) = v , c(0) = c . 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.

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 pri- mary 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 dif-fusion 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 predatorprey 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.
Acknowledgements Purnedu Mishra sincerely thanks the Faculty of Mathematics, Informatics and Mechanics, University of Warsaw Poland, for providing the research facilities. This work has been carried out during the tenure of an ERCIM, Alain Bensoussan Fellowship Programme. Sincere thanks to Prof. Dariusz Wrzosek and Prof. Malay Banerjee for helpful discussions during the revision process. Authors extend their appreciation to handling editor and all potential reviewers for their constructive comments.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.