Mathematics in Applied Sciences and Engineering https://ojs.lib.uwo.ca/mase Volume 4, Number 2, June 2023, pp.79-99 https://doi.org/10.5206/mase/15231 DIFFUSION-DRIVEN INSTABILITY AND PATTERN FORMATION IN A PREY-PREDATOR MODEL WITH FEAR AND ALLEE EFFECT DEBJIT PAL, DIPAK KESH, AND DEBASIS MUKHERJEE Abstract. This paper analyses a predator-prey model with Holling type II response function incorpo- rating Allee and fear effect in the prey. The model includes intra species competition among predators. We find out the local dynamics as well as Hopf bifurcation by considering level of fear as bifurcation parameter. The condition for diffusion-driven instability and patterns are then demonstrated in rela- tion to the system’s ecological parameters and diffusion coefficients. Intra-specific competition affects the dynamics of the system and Turing pattern formation. Moreover, output of results are verified through numerical simulation. Thus, from a dynamical standpoint, the considered model seems to be relevant in the field of ecology. 1. Introduction The investigation of diverse systems relating to the interaction between predators and preys, which are the building blocks of an ecosystem, is among the most intriguing research areas for ecologists, biologists, and even mathematicians. The most important aspect of a predator-prey interaction is that predator development is determined by its functional response to the prey species, which quantifies how much prey each predator consumes per unit of time. There are two techniques to determining a predator’s influence on prey in predator-prey interactions [3, 4]. One involves a predator killing its prey directly, while the other involves an indirect tactic such as fear of predation, Allee effect caused by predation, and so on. In comparison to direct predation, indirect predation has a more efficient and long-lasting evolutionary effect. In reality, all animals exhibit a variety of anti-predator responses, such as habitat use, alterations in hunting behaviour, physiology etc. In a predator-prey system, the well-established idea is that predators can only impact prey by destroying prey populations. Direct prey killing by predators is far more common in ecology compared to indirect influence. However, certain experimental findings demonstrate that, in addition to direct killing, fear of predator has a significant influence on prey behaviour and physiology [7]. Because of their fear of predators, preys are often forced to relocate their grazing region to a safer location by surrendering their most preferred regions. As a result, their chances of survival rise immediately, but their long-term fitness drops. The prey’s fecundity rate is reduced as a result of this. As an illustration of indirect dread Zanette et al. [36] investigated a field study on song sparrows over a full mating season, utilising electric fences and seines to deliberately prevent direct predation and managing predation risk through predator call playbacks. They studied the effect of fear on demography and discovered a 40% decrease in offspring generated each year in free-living song sparrows population. Sasmal and Takeuchi [17] investigated a predator- prey system including fear effect and explored the multi-stability of the system and the occurrence of Received by the editors 16 September 2022; accepted 13 March 2023; published online 2 April 2023. 2020 Mathematics Subject Classification. Primary 35B32, 35B36; Secondary 92D25. Key words and phrases. Allee effect, fear effect, diffusion, Hopf bifurcation, Turing instability, pattern. D.Pal was supported by CSIR (File No. 09/0096(12492)/2021-EMR-I), Government of India. 79 https://ojs.lib.uwo.ca/mase https://dx.doi.org/https://doi.org/10.5206/mase/15231 80 D. PAL, D. KESH, AND D. MUKHERJEE Hopf bifurcation. Xiao and Li [31] studied a predator-prey mutual interaction model with fear impact. Several work in predator-prey model was also done incorporating fear effect [16, 37, 40, 28, 20, 11, 33]. One of the most intriguing phenomena in ecology is the Allee effect, which was put forth by Allee [2]. The Allee effect is supported by the correlation between species density and relative growth rate at low population densities which have been observed in many natural populations, including plants, insects, marine invertebrates, birds, and mammals. Some of the factors underlying the Allee effect include difficulties in finding a mate, reproductive enablement, predation, inbreeding anxiety, environment conditioning, low density social disorder, and so on. Allee effects are broadly separated into two types: strong Allee effects and weak Allee effects. The strong Allee effect has a population density threshold that causes the species to become extinct if it falls below this population density [9]. The weak Allee effect, on the other hand, happens when the growth rate slows but remains positive at low population density. Ecologists have given close attention to this issue as it connects to species extinction. Sen and Banerjee [19] investigated several sorts of bifurcation behaviour in the Holling type II model with a strong Allee effect in growth function of prey and a density-dependent mortality rate for predators. With the use of the Holling type-III functional response, Vishwakarma and Sen [25] looked into the influence of the Allee effect on prey and hunting cooperation in predators. Several scientists have also examined the Allee effect in predator-prey models [1, 18, 35]. Many scholars have shown interest in predator-prey model incorporating both fear effect and Allee effect. Liu [13] considered a predator-prey model that included fear and Allee effects and studied qualitative nature of the system. Lie et al. [12] explored the global stability, saddle-node bifurcation, transcritical bifurcation, and Hopf bifurcation in a predator-prey model, taking into account the fear effect and additive Allee effect. Xie [32] discussed the influence of the Allee and fear effects on a Holling type-II predator-prey model by considering the system as  du dt = u ( ru (a + u) (1 + fv) −d1 −ku ) − euv 1 + bu , dv dt = −d2v + ceuv 1 + bu . (1.1) where u (t) and v (t) are the respective densities of prey and predator populations at time t. As Allee effect also acts on prey’s birth rate so birth rate of prey is expressed as Q(a,u) = ru a+u , where u denotes density of prey at time t (u > 0), r(> 0) denotes prey’s maximum birth rate and a(> 0) presents the level of Allee and Q(a,u) satisfied limu→0 Q(a,u) = 0, limu→+∞Q(a,u) = r, lima→0 Q(a,u) = r, lima→+∞Q(a,u) = 0 and ∂Q(a,u) ∂a < 0. Since certain experimental tests revealed that the fear effect can significantly slow down the repro- duction process so prey’s birth rate, Q (a,u) is scaling by a term N(f,v) = 1 (1+fv) , this speaks of the price of an anti-predator response brought on by fear, and the parameter f indicates fear level. Con- sidering the experimental results and biological viewpoint, the factor N(f,v) must satisfy the following conditions: (1) N (0,v) = 1 : there will be no reduction in prey production if prey exhibits no anti-predator activities(fear of predation). (2) N (f, 0) = 1 : there is no decrease in prey production due to anti-predator behaviours if there is no predator. (3) limf→+∞N(f,v) = 0 : prey production decreases to zero when anti-predator behaviour is extremely high. (4) limv→+∞N(f,v) = 0 : due to intense anti-predator responses, if the predator population is sufficiently large, prey reproduction is decreased to zero. (5) ∂N(f,v) ∂f < 0 : the production of prey declines as anti-predator behaviours rise. DIFFUSION-DRIVEN INSTABILITY AND PATTERN FORMATION IN A PREY-PREDATOR SYSTEM 81 (6) ∂N(f,v) ∂v < 0 : prey productivity declines as the size of the predator population rises. Inspired by the preceding model (1.1), in a prey-predator system with a Holling type-II functional response, we wish to examine the combined impact of include both fear and the Allee effect in addition to predator’s intra-specific competition which is added with model (1.1). The new model can be written as ;   du dt = u ( ru (a + u) (1 + fv) −d1 −ku ) − euv 1 + bu , dv dt = −d2v −hv2 + ceuv 1 + bu . (1.2) where it is assumed that all of the parameters are positive and are shown in table 1. Parameters Meaning value u prey density at time t − v predator density at time t − a level of Allee 0.1 b handling time of predators 0.9 c conversion efficiency 0.6 d1 mortality of prey 0.9 d2 mortality of predator 0.2 e consumption rate of predator 0.85 k intraspecific competition among the prey species 0.1 h intraspecific competition among the predator species − r maximum birth rate of prey − f level of fear induced by predator − Table 1. Description of parameters and values of some parameters in model (1.1) and (1.2) . The idea of spatial pattern development in predator-prey systems by Turing instability plays a crucial role in ecology as well as other real-world issues in science, engineering, and technology. Turing insta- bility refers to the fact that diffusion converts a stable system to an unstable one. Patterns are formed by discontinuities in population distribution over time and space, and they reveal environmental irreg- ularities. Diffusion significantly affects the development of patterns. It is self-evident that spatial and spatiotemporal pattern formation is critical in ecology. The patterns we see in biological and ecological communities are those that are required for their survival. However, the processes that produce patterns are still unknown. The reaction-diffusion hypothesis of pattern creation was initially put forward by Alan Turing in 1952 in his seminal study on the chemical basis of morphogenesis [24]. Later, Wolpert [29] suggested a phenomenological idea of pattern creation and gave a very concise, non-technical ex- planation of how animal pattern generation occurs [30], which has resulted in an immense number of enlightening experimental research, many of which are also related to marine ecology. Chen and Wang [6] investigate qualitative features of a predator-prey model with diffusion system and demonstrate the presence and absence of positive equilibrium point. In recent times, many researchers give attention in the field of pattern formation [15, 5, 21, 34, 10, 23, 27] such as Sun et al. have studied several dynamics of pattern in a spatial prey predator model incorporating Allee [22], Zhang et.al. examined diffusive ratio dependent predator-prey model taking delay in fear effect and discussed the different bifurcations 82 D. PAL, D. KESH, AND D. MUKHERJEE and patterns induced by delay [38], Duan et.al. [8] looked into Hopf-Hopf bifurcation and chaotic at- tractors considering fear effect in delayed diffusive predator-prey model. Zhang, Zhao and Yuan [39] explored the effect of discontinuous harvesting in diffussive pre-predator model in presence of Allee and fear effect. Jun Ma et al. [14, 26] investigate an interesting study of patterns in neuronal network. In general, reaction-diffusion equations are frequently used to represent spatio-temporal dynamics, and the study of such systems adds to our understanding of the performance of spatio-temporal instabilities in ecosystems with a large degrees of freedom. When we extend model (1.2) by integrating random dispersal of both species, the diffusive model becomes:   ∂u ∂t = u ( ru (a + u) (1 + fv) −d1 −ku ) − euv 1 + bu + D1∇2u, ∂v ∂t = −d2v −hv2 + ceuv 1 + bu + D2∇2v. (1.3) Considering the two- dimensional spatial domain Ω = [0,L] × [0,L] ∈ R2, ∇2 = ∂ 2 ∂x2 + ∂ 2 ∂y2 and subject to non-negative initial conditions and zero flux boundary conditions. Here for ecological feasibility, prey and predator’s respective diffusion coefficients D1 and D2 are choosen to be positive, and u (x,y,t) and v (x,y,t) denote the respective densities of prey and predator, at spatial location (x,y) ∈ Ω and time instant t. All the other parameters of the model are described in previous section in table 1. In this paper, for simplicity of calculation, birth rate of the prey, r is considered to be constant, not density dependent. The prime objective of our analysis is to examine the continuous model without diffusion as well as diffusion-driven instability and pattern creation in 2D coordinates. This article is arranged as follows: In section 2, we address the existence and local stability of the equilibrium points. The occurrence of Hopf bifurcation around the interior equilibrium point is established in section 4. Section 5 reveals numerical simulations of the model in absence of diffusion. In section 6, we explore the model with diffusion with initial and boundary conditions and in two subsections we establish the condition for diffusion-driven instability and the numerical simulation of pattern formation. Finally, the article is ended in section 7 with a discussion where a short comparison between the presence and absence of intra-specific competition among predator is also shown. 2. Equilibrium Points Here we shall find all the equilibrium points of system (1.2). Firstly, E0 (0, 0) is always an equilibrium point of system (1.2) . Searching for predator free equilibrium points, let us take (u′, 0) be such an equilibrium point. Then, u′ is the positive root of ku2 + (d1 + ak −r) u + ad1 = 0. The solution is, u = −(d1 + ak −r) ∓ √ ∆1 2k . (2.1) where, ∆1 = (d1 + ak −r) 2 − 4akd1. When (d1 + ak −r) < 0 and ∆1 > 0 i.e. 0 < a < a1 = (√ r − √ d1 )2 k DIFFUSION-DRIVEN INSTABILITY AND PATTERN FORMATION IN A PREY-PREDATOR SYSTEM 83 holds, then model (1.2) has two predator free or axial equilibrium points i.e. two values of u′, given by E1 (u1, 0) = ( −(d1 + ak −r) − √ ∆1 2k , 0 ) , E2 (u2, 0) = ( −(d1 + ak −r) + √ ∆1 2k , 0 ) . An interior equilibrium point E∗ (u∗,v∗) will be the positive solution of the two equations given by:  ru (a + u) (1 + fv) −d1 −ku− ev 1 + bu = 0, −d2 −hv + ceu 1 + bu = 0. (2.2) From second equation of (2.2), v = Mu−d2 h(1 + bu) . (2.3) where M = ce− bd2. To find u, putting the expression of v from (2.3) into the first equation of (2.2) we get a polynomial of degree 5 in u, in the form (EG) u5 + ( EF + DG−rb3h2 ) u4 + ( CG + DF − 3rb2h2 ) u3 + ( BG + CF − 3rbh2 ) u2. + ( AG + BF −rh2 ) u + AF = 0. (2.4) where, A = a (d1h−ed2) B = 2abhd1 + ahk + hd1 + Mea−ed2, C = 2abhk + ab2hd1 + 2bhd1 + hk + Me, D = b2hd1 + ab 2hk + 2bhk, E = b2hk, F = h−fd2, G = bh + fM. We know that a polynomial of degree 5 must have a real root, positive or negative. If constant term AF < 0 i.e. either A > 0 and F < 0 or A < 0 and F > 0 ⇒ either ed2 d1 < h < fd2 or ed2 d1 > h > fd2 , then the equation (2.4) must have at least one positive root, say u∗. Putting that u∗ in Equation (2.3) we can get positive v = v∗(say) if Mu∗ − d2 > 0 ⇒ u∗ > d2M and M > 0 ⇒ ce > bd2 and an interior equilibrium point of given model (1.2) exists only if both u∗ and v∗ are positive. Again one can also give conclusive evidence about the existence of the interior equilibrium point from graphs of isoclines as shown in Fig 1. 3. Local Stability Analysis of Equilibrium Points: We first calculate the Jacobian matrix of system (1.2) at any arbitrary point (u,v) and it is given by; J(u,v) =   ru(2a+u) (a+u)2(1+fv) −d1 − 2ku− ev(1+bu)2 −rfu2 (a+u)(1+fv)2 − eu (1+bu) cev (1+bu)2 −d2 − 2hv + ceu(1+bu)   . (3.1) Now we shall discuss the local stability of the system (1.2) near each of the equilibrium point. For E0 (0, 0), from (3.1) we get 84 D. PAL, D. KESH, AND D. MUKHERJEE Figure 1. Position of x-isoclines(red) and y-isoclines(blue) and interior equilibrium point is indicated by a circle, for a particular set of values of parameters. J (0, 0) =   −d1 0 0 −d2   . The eigenvalues are −d1 and −d2, both are real negative hence extinction steady state E0 (0, 0) become stable i.e., when concentration of prey and predator is belonged in the attraction region of E0, both the populations will die out. Now from (3.1) for E1 (u1, 0) we get, J(u1, 0) =   ru1(2a+u1) (a+u1) 2 −d1 − 2ku1 −rfu21 (a+u1) − eu1 (1+bu1) 0 −d2 + ceu1(1+bu1)   . (3.2) Now from equation (2.1) we can write√ ∆1 = r −d1 −ak − 2ku1, =⇒ u1 √ ∆1 = ad 2 1 −ku 2 1. Using ru ′ a+u′ −d1 −ku′ = 0, and using the last expression, the Jacobian takes the form J(u1, 0) =   u1 √ ∆1 (a+u1) −rfu21 (a+u1) − eu1 (1+bu1) 0 −d2 + ceu1(1+bu1)   . (3.3) The two eigenvalues are given by λ1 = u1 √ ∆1 (a + u1) > 0, λ2 = −d2 + ceu1 (1 + bu1) . DIFFUSION-DRIVEN INSTABILITY AND PATTERN FORMATION IN A PREY-PREDATOR SYSTEM 85 So λ2 > 0 if −d2 + ceu1(1+bu1) > 0 ⇒ b < b1 = ( ce d2 − 1 u1 ) ; and λ2 < 0 if b > b1. Thus, E1 (u1, 0) is a unstable node for 0 < b < b1 and E1 (u1, 0) is a saddle point for b > b1. Similarly for E2 (u2, 0), we get, J(u2, 0) =   −u2 √ ∆1 (a+u2) −rfu22 (a+u1) − eu2 (1+bu2) 0 −d2 + ceu2(1+bu2)   . (3.4) The two eigenvalues are given by λ1 = −u2 √ ∆1 (a + u2) < 0, λ2 = −d2 + ceu2 (1 + bu2) . So λ2 > 0 if −d2 + ceu2(1+bu2) > 0 =⇒ b < b2 = ( ce d2 − 1 u2 ) , and λ2 < 0 if b > b2. Therefore, E2 (u2, 0) is a saddle point for 0 < b < b2 and E2 (u2, 0) is a stable node for b > b2. Summarizing the above, for the predator free equilibrium points E1 and E2 , the results are repre- sented in table2 Existence criterion Equilibria Nature 0 < a < a1 E1 (u1, 0) unstable node if 0 < b < b1 and saddle if b > b1 0 < a < a1 E2 (u2, 0) saddle if 0 < b < b2 and locally asymp- totically stable if b > b2 Table 2. Existence and stability condition of predator extinction steady state of model (1.2). . From the above discussion we can conclude that if Allee effect is very weak i.e. value of the parameter a is smaller than threshold value a1, then predator free equilibrium points exit. In addition to it also if the handling time of predator b is rather greater than a threshold value, then E1 is a saddle point (b > b1) and E2 is locally asymptotically stable (b > b2). The quadrant I is separated into two regions, attraction regions of E0 and E2, respectively. On the other hand if the handling time of predator b is rather smaller than a threshold value, then E1 is unstable (0 < b < b1) and E2 is saddle (0 < b < b2). If system (2) has no coexisting equilibrium point then based on on the level of Allee, handling time of the predator and initial population density, either predator species will become extinct (when trajectory tends to E2 ) or both prey and predator populations will go to extinct (when trajectory tends to E0). Now to determine the stability of E∗ (u∗,v∗), we have to consider the sign of determinant and trace of J at (u∗,v∗) . It is known that the system will be stable if tr < 0 and det > 0. From (3.1) for E∗ (u∗,v∗), we get, J (u∗,v∗) =   ru∗(2a+u∗) (a+u∗)2(1+fv∗) −d1 − 2ku∗ − ev ∗ (1+bu∗)2 −rfu∗2 (a+u∗)(1+fv∗)2 − eu ∗ (1+bu∗) cev∗ (1+bu∗)2 −d2 − 2hv∗ + ceu ∗ (1+bu∗)   =:   J11 J12 J21 J22   (3.5) 86 D. PAL, D. KESH, AND D. MUKHERJEE where J12 = −rfu∗2 (a + u∗) (1 + fv∗) 2 − eu∗ (1 + bu∗) < 0, J21 = cev∗ (1 + bu∗) 2 > 0, J22 = −d2 − 2hv∗ + ceu∗ (1 + bu∗) = −hv∗ < 0 (By second equation of (2.2) ), J11 = ru∗ (2a + u∗) (a + u∗) 2 (1 + fv∗) −d1 − 2ku∗ − ev∗ (1 + bu∗) 2 = u∗ ( ra (a + u∗) 2 (1 + fv∗) −k + bev∗ (1 + bu∗) 2 ) (By first equation of (2.2) ). The characteristic equation is, λ2 − (J11 + J22) λ + (J11J22 −J12J21) = 0. where det J= (J11J22 −J12J21) and trJ = (J11 + J22) = v ∗ ( ae (1 + bu∗) (a + u∗) + beu∗ (1 + bu∗) 2 ) − k(u∗)2 −ad1 a + u∗ −hv∗ by equation (2.2). We first consider the stability condition on trJ. For stability of the system we must have trJ < 0, that is, v∗ ( ae (1 + bu∗) (a + u∗) + beu∗ (1 + bu∗) 2 ) < ( k(u∗)2 −ad1 a + u∗ + hv∗ ) =⇒ v∗ ( ae (1 + bu∗) (a + u∗) + beu∗ (1 + bu∗) 2 ) + ad1 a + u∗ < ( k(u∗)2 a + u∗ + hv∗ ) . (3.6) Then from the equation (3.6), we have e < e∗ where e∗ = (1 + bu∗) 2 [ (a + u∗) hv∗ + ( ku∗2 −ad1 )] v∗ (a + 2abu∗ + bu∗2) . (3.7) Thus, trJ < 0 if e < e∗. We next consider the stability condition on detJ. It is known that for E∗ to be stable, we must have detJ > 0, that is, ⇒ J11J22 −J12J21 > 0. (3.8) Since (–J12J21) > 0 , so equation (3.8) holds good in two cases below. Case-I: J11J22 > 0. Since J22 < 0, the above holds if J11 < 0, that is, rau∗ (a + u∗) 2 (1 + fv∗) + beu∗v∗ (1 + bu∗) 2 < ku∗. (3.9) Case-II: J11J22 < 0 and in addition |J11J22| < |J12J21|. Since J22 < 0, the above holds if J11 > 0 and |J11J22| < |J12J21|. =⇒ rau ∗ (a+u∗)2(1+fv∗) + beu ∗v∗ (1+bu∗)2 > ku∗ and |J11J22| < |J12J21|. (3.10) DIFFUSION-DRIVEN INSTABILITY AND PATTERN FORMATION IN A PREY-PREDATOR SYSTEM 87 Summarizing the above conditions on trJ and det J, we conclude that E∗ (u∗,v∗) is stable when e < e∗ and equation (3.9) or (3.10) hold. 4. Hopf Bifurcation In the study of continuous prey-predator models, bifurcation is possibly the most crucial topic. Here, we use the predator-induced fear level, f, as a bifurcation parameter in the investigation of the Hopf bifurcation of the system (1.2). For Hopf bifurcation, we need to find out a critical value of f = f1 such that (1) trJ = J11 + J22 = 0 at f = f1. (2) detJ = J11J22–J12J21 > 0 at f = f1. The characteristic equation of J (u∗,v∗) takes the form λ2 + detJ = 0 at f = f1 whenever E∗ (u∗,v∗) is feasible. As a result, two eigenvalues are conjugate to one another and are purely imaginary. (3) If the eigenvalues of (3.1) around E∗ (u∗,v∗) have the general form λ = λ1 +iλ2 , then dλ1 df |f=f1 6= 0. The characteristic equation at E∗ is given by λ2 − (trJ) λ + detJ = 0. Equating the real and imaginary parts of the both sides, we get( λ21 −λ 2 2 ) − (trJ) λ1 + detJ = 0, 2λ1λ2 − (trJ) λ2 = 0. (4.1) Differentiating the aforementioned equation with respect to the bifurcation parameter f, we obtain (2λ1 − trJ) dλ1 df − 2λ2 dλ2 df −λ1 d df (trJ) + d df (detJ) = 0, 2λ2 dλ1 df + (2λ1 − trJ) dλ2 df −λ2 d df (trJ) = 0. (4.2) Now eliminating dλ2 df from the above two equation of (4.2) we obtain dλ1 df as dλ1 df = 1 P 2 + 4λ22 [ P ( λ1 d df (trJ) − d df (detJ) ) + 2λ22 d df (trJ) ] , where P = (2λ1 − trJ). Using the facts trJ|f=f1 = 0 and λ1|f=f1 = 1 2 trJ|f=f1 = 0, we have, dλ1 df |f=f1 = 1 2 d df (trJ) |f=f1 = 1 2 d df [( u∗ ( ra (1 + fv∗) (a + u∗) 2 + bev∗ (1 + bu∗) 2 ) −k ) −hv∗ ] f=f1 6= 0. As a result, the system admits a Hopf bifurcation near E∗ (u∗,v∗) at f = f1, where f1 is the positive solution of the equation[( u∗ ( ra (1 + fv∗) (a + u∗) 2 + bev∗ (1 + bu∗) 2 ) −k ) −hv∗ ] f=f1 = 0. In this section Jacobian J is calculated at E∗ (u∗,v∗). 88 D. PAL, D. KESH, AND D. MUKHERJEE 5. Numerical Simulation Now we examine the analytical results through numerical simulation by taking different values of parameters. Here we fix seven parameters out of ten as b = 0.9, c = 0.6, e = 0.85, a = 0.1, k = 0.1, d1 = 0.9, d2 = 0.2. Example 1 : Set h = 0.18, f = 1.2, r = 1. For these choice of parameters E0 (0, 0) becomes stable (Fig. 2). (a) Phase diagram of prey and predator (b) Phase diagram of prey and predator Figure 2. For h = 0.18, f = 1.2, r = 1 with initial value (1, 0.3) Example 2 : Choose h = 0.18, f = 1.2, r = 1.5. Then 0 < a < a1 = 2.7606 which implies that two predator free equilibrium points E1 (0.1567, 0) and E2 (5.7433, 0) exist along with the extinct equilibrium point E0 (0, 0). Since b1 = −3.8316 and b2 = 2.3759 so from table 2 we conclude that E1 and E2 both are saddle points hence the trajectory goes to (0, 0) (Fig. 3). Moreover if we increase the value of k from 0.1 to 0.7 then 0 < a < a1 = 0.3944. Hence two predator free equilibrium points E1 (0.2571, 0) and E2 (0.5, 0) exist along with the extinct equilibrium point E0 (0, 0). Here b1 = −1.3395 and b2 = 0.55 so table 2 concludes that E0 and E2 are locally asymptotically stable and E1 is saddle. The attraction zones of E0 and E2 divide the first quadrant into two sections. Both populations will become extinct as the initial population of prey lies in the attraction region of E0. Whereas the predator population will only begin to decline when initial population of prey lies in the attraction region of E2 instead of E0. DIFFUSION-DRIVEN INSTABILITY AND PATTERN FORMATION IN A PREY-PREDATOR SYSTEM 89 (a) Phase diagram of prey and predator (b) Phase diagram of prey and predator Figure 3. For h = 0.18, f = 1.2, r = 1.5, k = 0.1 with initial value (1, 0.01) (a) Phase diagram of prey and predator (b) Phase diagram of prey and predator Figure 4. For h = 0.18, f = 1.2, r = 1.5, k = 0.7 with initial value (1, 0.01) (blue) and (1, 0.1) (red). Now if we also fixed h = 0.2, r = 2.2 then for increasing values of level of fear f, the following dynamics will obtain. Example 3 : For f = 1.05 the system (1.2) has an equilibrium point E∗(1.4339, 0.5963). Here e∗ = 1.1893, hence e < e∗ so E∗ becomes stable (Fig. 5). 90 D. PAL, D. KESH, AND D. MUKHERJEE (a) Phase diagram of prey and predator (b) Phase diagram of prey and predator Figure 5. For h = 0.2, f = 1.05, r = 2.2 with initial value (1, 0.4) Example 4: For f = 1.84 the system (1.2) possesses unique equilibrium point E∗(1.0658, 0.3872). In this case, the criterion of Hopf bifurcation are fulfilled i.e. trJ = 0 at f = f1, detJ > 0 at f = f1 and dλ1 df |f=f1 6= 0 (f1 = 1.84) and the model admits stable limit cycle with increasing f (Figs. 6 and 7). (a) Phase diagram of prey and predator (b) Phase diagram of prey and predator Figure 6. For h = 0.2, f = 1.84, r = 2.2 with initial value (1, 0.4) DIFFUSION-DRIVEN INSTABILITY AND PATTERN FORMATION IN A PREY-PREDATOR SYSTEM 91 (a) prey(u) vs fear effect(f) (b) predator(v) vs fear effect(f) (c) 3-D picture Figure 7. Hopf bifurcation Example 5 : For f = 2.2 the system (1.2) has unique equilibrium point E∗(0.9906, 0.3354). Here e∗ = 0.7466 hence e > e∗ and E∗ is an unstable source and the trajectory goes to a limit cycle (Fig. 8). (a) Phase diagram of prey and predator (b) Phase diagram of prey and predator Figure 8. For h = 0.2, f = 2.2, r = 2.2 with initial value (1, 0.365) 92 D. PAL, D. KESH, AND D. MUKHERJEE 6. Model With Diffusion In this part, we consider the model (1.3) and examine the impact of diffusion about E∗ (u∗,v∗), subject to non-negative initial conditions: u (x,y, 0) = u0 (x,y) and v (x,y, 0) = v0 (x,y) for (x,y) ∈ Ω, and zero flux boundary conditions:( ∂ ∂x u (x,y,t) ) |x=0 = 0, ( ∂ ∂x u (x,y,t) ) |x=L = 0,( ∂ ∂y u (x,y,t) ) |y=0 = 0, ( ∂ ∂y u (x,y,t) ) |y=L = 0,( ∂ ∂x v (x,y,t) ) |x=0 = 0, ( ∂ ∂x v (x,y,t) ) |x=L = 0,( ∂ ∂y v (x,y,t) ) |y=0 = 0, ( ∂ ∂y v (x,y,t) ) |y=L = 0. 6.1. Diffusion-driven Instability. The formation of patterns can occur for a variety of reasons, but the focus of this article is on the pattern formation caused by diffusion-driven instability. The criterion for stability of uniform steady state E∗ (u∗,v∗) without diffusion are given below: tr (J (u∗,v∗)) = (J11 + J22) < 0 and det (J (u ∗,v∗)) = J11J22–J12J21 > 0. Now to determine the Turing instability we linearize the system (1.3) ( We expand u and v in spatial domain as u (x,y,t) = u∗ + u1 (x,y,t) and v (x,y,t) = v ∗ + v1 (x,y,t) where |u1| << u∗ and |v1| << v∗ and ignoring all the non linear terms in u and v). Patterns are time independent and spatially heterogeneous solution of system (1.3). Let us take the solution is separable, so take u1 (x,y,t) = A (t) e i(kxx+kyy), and v1 (x,y,t) = B (t) e i(kxx+kyy), where k = √ k2x + k 2 y denotes the wave number of the solution. The additional two diffusive terms which are incorporated in model (1.2) to get model (1.3) are D1 ( ∂2 ∂x2 ( A (t) ei(kxx+kyy) ) + ∂2 ∂y2 ( A (t) ei(kxx+kyy) )) = −D1k2 ( A (t) ei(kxx+kyy) ) and D2 ( ∂2 ∂x2 ( B (t) ei(kxx+kyy) ) + ∂2 ∂y2 ( B (t) ei(kxx+kyy) )) = −D2k2 ( B (t) ei(kxx+kyy) ) . Then the corresponding Jacobian matrix for the system (1.3) is given by J (k) =   J11 −D1k2 J12 J21 J22 −D2k2   . The characteristic equation is, µ2 − tr (J (k)) µ + det (J (k)) = 0. where tr (J (k)) = J11 −D1k2 + J22 −D2k2 = tr (J (u∗,v∗)) − (D1 + D2) k2 DIFFUSION-DRIVEN INSTABILITY AND PATTERN FORMATION IN A PREY-PREDATOR SYSTEM 93 and det (J (k)) = ( J11 −D1k2 )( J22 −D2k2 ) −J12J21 = det (J (u∗,v∗)) − (D1J22 + D2J11) k2 + D1D2k4. Applying Routh-Hurwitz criteria from the characteristic equation we can conclude that E∗ (u∗,v∗) is locally asymptotically stable for all k > 0 if tr (J (k)) = tr (J (u∗,v∗)) − (D1 + D2) k2 < 0 and det (J (k)) = det (J (u∗,v∗)) − (D1J22 + D2J11) k2 + D1D2k4 > 0. From the stability of diffusion free system we get that tr (J (u∗,v∗)) < 0 and since D1,D2 and k are all positive so tr (J (k)) < 0 always. Since in case of forming Turing pattern the system have to be unstable with respect to spatially heterogeneous perturbations. So the condition to induce Turing instability is that, for some k > 0, det (J (k)) = det (J (u∗,v∗)) − (D1J22 + D2J11) k2 + D1D2k4 < 0. (6.1) Now we take derivative to det (J (k)) with respect to k2 and equate to zero so that we can obtain following critical value k2c = D1J22 + D2J11 2D1D2 . To get such feasible critical value k2c , the necessary condition is D1J22 + D2J11 > 0. For this critical value k2c , det (J (k)) will be minimum and that minimum value is obtained by putting the value of k2c in the expression of det (J (k)) in (6.1) and that minimum value is thus obtained is, det (J (kc)) = det (J (u ∗,v∗)) − (D1J22 + D2J11) 2 4D1D2 . For Turing instability the minimum value of det (J (k)) i.e det (J (kc)) must also be negative when D1J22 + D2J11 > 2 √ D1D2 (det (J (u∗,v∗))) =⇒ D1J22 + D2J11 > 2 √ D1D2 (J11J22–J12J21). Now we are in a position to list the conditions for diffusion-driven instability and the list is as follows: (1) J11 + J22 < 0, (2) J11J22–J12J21 > 0, (3) D1J22 + D2J11 > 0, (4) D1J22 + D2J11 > 2 √ D1D2 (J11J22–J12J21). 6.2. Numerical simulation. In this section, in support of our analytic and theoretical finding, we describe numerical observations of system (1.3) in Ω ∈ R2. To simulate system (1.3), we must employ the Forward Time Centered Space (FTCS) technique to convert an infinite dimensional continuous model to a finite dimensional discrete model. Simulation is performed to the model (1.3) for all (x,y) ∈ Ω = [0,L] × [0,L] where δt = 0.05 (time step size) and δx = 0.5, δy = 0.5 (space step sizes) and initial data with a random perturbations around E∗ (u∗,v∗). Numerical trials are conducted till the pattern achieves a stationary state, at which time the system’s behaviour ceases to change. By studying several sorts of spatial dynamics using numerical evaluations, we discovered that the distributions of both species are consistently identical. Other parameters value are taken as d1 = 0.9, d2 = 0.2, k = 0.1, b = 0.9, e = 0.85, c = 0.6, r = 2.2, a = 0.1, h = 0.2, D1 = 0.01, D2 = 1.2 and different input parameter f are chosen as 1.05, 1.20, 1.35, 1.50, 1.65 and 1.80 and in all the cases the conditions of diffusion driven instability are satisfied. Figs. 9 and 10 display respective spatial patterns of prey and predator species around E∗ (u∗,v∗) in 2D space. An interesting finding from Figs. 9 and 10 is that the bounds for prey 94 D. PAL, D. KESH, AND D. MUKHERJEE population u and predator population v over 2D space is varying with the change of system parameter f. Spots pattern (H -pattern) is generated for prey species, while for predator species, larger values of f (f = 1.65,f = 1.8) form a spot dominated combination of spots and very short stripes. We now provide some biological justifications for the trend shown in the aforementioned figures. For prey, the spots pattern consists of blue (low population density) quadrilaterals on a red (high population density) backdrop, indicating that the prey species are concentrated in isolated areas having low population density, whereas the rest of the region is densely populated. (a) f = 1.05 (b) f = 1.2 (c) f = 1.35 (d) f = 1.5 (e) f = 1.65 (f) f = 1.8 Figure 9. Evolution of prey species for d1 = 0.9, d2 = 0.2, k = 0.1, b = 0.9, e = 0.85, c = 0.6, r = 2.2, a = 0.1, h = 0.2, D1 = 0.01, D2 = 1.2 and different values of f. 7. Discussion The present investigation deals with a predator-prey system with Allee and fear effects on prey and Holling type-II functional response. Intra-specific competition among predators is also taken under consideration. We extend the concept of generalized predator–prey model to spatially extended systems, in which the motion can be described by a diffusion term to the system (1.2). Firstly, we have studied the model without diffusion and have found the equilibrium points and deduce that interior equlibrium point is stable when consumption rate of predator, e less than a threshold value, e∗ and equation (3.9) or (3.10) hold. Then we have examined the existence of hopf bifurcation and limit cycle taking level of fear, f as bifurcation parameter (Figs. 5- 8). Xie [32] discussed system (1.1) but the difference between the system (1.1) and (1.2) is that intra- specific competition among predators, h is considered in system (1.2). Consideration of h does not affect the trivial and axial equilibrium points, it only affects the stability of interior equilibrium point E∗ and it is illustrated by an example. If we take h = 0 and the other parameters value same as example 3 in section 5, then simulation result, Fig. 11 shows that if h is not taken into consideration (for system DIFFUSION-DRIVEN INSTABILITY AND PATTERN FORMATION IN A PREY-PREDATOR SYSTEM 95 (a) f = 1.05 (b) f = 1.2 (c) f = 1.35 (d) f = 1.5 (e) f = 1.65 (f) f = 1.8 Figure 10. Evolution of predator species for d1 = 0.9, d2 = 0.2, k = 0.1, b = 0.9, e = 0.85, c = 0.6, r = 2.2, a = 0.1, h = 0.2, D1 = 0.01, D2 = 1.2 and different values of f. (a) For h = 0.2 (b) For h = 0 Figure 11. The effect of intra-specific competition among predators on the stability of interior equilibrium point. For f = 1.05, r = 2.2, b = 0.9, c = 0.6, e = 0.85, a = 0.1, k = 0.1, d1 = 0.9, d2 = 0.2. (1.1) ), E∗ becomes unstable and trajectory goes to extinct. But when h is considered (for system (1.2)), for the same value of parameters, E∗ becomes stable and both species exist simultaneously. So the stability of E∗ is affected when intra-specific competition among predators is considered. Another numerical simulation in Fig. 12 shows that for various initial populations, in absence of h the trajectory 96 D. PAL, D. KESH, AND D. MUKHERJEE always goes to (0, 0) i.e. both the population goes to extinction but in presence of h, for some initial populations (red coloured region in Fig. 13) trajectory goes to (0, 0) and for some initial populations ( blue coloured region in Fig. 13) trajectory goes to E∗(1.4339, 0.5963) and both the populations survive. So the region is divided into two parts, attraction region of E0 and attraction region of E ∗. Hence interior equilibrium point becomes stable when intra-specific competition among predators is considered. (a) For h = 0.2 (b) For h = 0 Figure 12. Role of h on stability of interior equilibrium point for the parameters set in example 3 of section 4 with various initial points (1, 0.4), (0.8, 0.8), (0.3, 0.8), (1, 1) and (1.9, 0.9). We then introduced diffusion in the model and spatial pattern generation mentioned in Sect. 5, the conditions necessary for diffusion-driven instability and stimulating Turing pattern formation cor- responding to the reaction-diffusion system have been achieved successfully. We have explored the presence of Turing patterns by varying the control parameter f. The contour image of the spatial patterns for preys and predators are depicted in Figs. 9 and 10 which show the temporal evolution of interacting populations at t = 10000 for various values of f. The color-bars represent that the den- sity of the prey population drops with rising fear level which supports the fact that as level of fear increases, prey moving away from the predator, but all prey responds to predation danger in different ways, including habitat changes, feeding, reduced mating, vigilance, and various physiological changes. Such anti-predator efforts may lower reproduction in the long run, but they are immediately helpful since they raise the likelihood of adult survival. In addition, fearful prey tend to forage less, which reduces their growth and reproduction rate of prey. As a result due to less-availability of food, predator population density also decreases. We observe that intra-specific competition among predators has also crucial role in the formation of Turing patterns. Conditions (1) and (3) of section 6.1 indicate that J11 and J22 must be of opposite signs for Turing instability. But in absence of h that is for system (1.1), J22 = 0 hence diffusion driven instability does not occur in system (1.1) in presence of self-diffusion only. But in presence of h Turing patterns appear for suitable choice of parameters. Hence h operates as an underlying mechanism capable of generating non-uniform spatial distributions of predators and prey via diffusion-driven instability. But the considered model (1.2) also have some limitations which also give us scopes for future works such as birth rate of the prey, r is considered to be constant, not density dependent. These kind of DIFFUSION-DRIVEN INSTABILITY AND PATTERN FORMATION IN A PREY-PREDATOR SYSTEM 97 Figure 13. Basin of attraction of equilibrium points E0 (Red coloured region) and E∗ (Blue coloured region) for parameter values b = 0.9, k = 0.1, r = 2.2, a = 0.1, e = 0.85, h = 0.2, f = 1.05, d1 = 0.9, d2 = 0.2, c = 0.6. assumptions are proper for some unicellular and bacterial species. The species which suffer from the fear effect and the population suffers from decline in their growth, growth rate is density dependent. So one can look into the dynamics of the system considering density dependent growth rate of prey. In the non-spatial model, we use 1 1+fv to represent the fear effect, however, in a diffusive model, consid- ering this term is not completely appropriate as the population is always moving in space, therefore, by considering that species are heterogeneous in space, predator density v changes so level of fear f is logically density dependent. But for sake of simplicity, we have taken level of fear as constant, neglecting the density dependence of f. So one can study the model taking density dependent fear. As a future work, to accurately evaluate the stability of various types of Turing patterns, the amplitude equations of Turing patterns may be computed. Acknowledgement. The authors are grateful to the anonymous reviewers and Associate Editor for their helpful comments for improving this paper. References 1. S. M. Al-Momen and R. K. Naji, The dynamics of sokol-howell prey-predator model involving strong allee effect, Iraqi Journal of Science (2021), 3114–3127. 2. W. C. Allee, Co-operation among animals, American Journal of Sociology 37 (1931), no. 3, 386–398. 3. J. Bhattacharyya and S. Pal, Stage-structured cannibalism with delay in maturation and harvesting of an adult predator, Journal of Biological Physics 39 (2013), no. 1, 37–65. 4. J. Bhattacharyya and S Pal, Hysteresis in coral reefs under macroalgal toxicity and overfishing, Journal of Biological Physics 41 (2015), 151–172. 5. Y. Cai, Z. Gui, X. Zhang, H. Shi, and W. Wang, Bifurcations and pattern formation in a predator–prey model, International Journal of Bifurcation and Chaos 28 (2018), no. 11, 1850140. 6. B. Chen and M. Wang, Qualitative analysis for a diffusive predator–prey model, Computers & Mathematics with Applications 55 (2008), no. 3, 339–355. 7. W. Cresswell, Predation in bird populations, Journal of Ornithology 152 (2011), no. 1, 251–263. 8. D. Duan, B. Niu, and J. Wei, Hopf-hopf bifurcation and chaotic attractors in a delayed diffusive predator-prey model with fear effect, Chaos, Solitons & Fractals 123 (2019), 206–216. 98 D. PAL, D. KESH, AND D. MUKHERJEE 9. S. N. Elaydi and R. J Sacker, Population models with allee effect: A new model, Journal of Biological Dynamics 4 (2010), no. 4, 397–408. 10. D. Geng, H. Wang, and W. Jiang, Nonlocal competition and spatial multi-peak periodic pattern formation in diffusive holling-tanner predator-prey model, Journal of Dynamics and Differential Equations (2022), 1–30. 11. R. P. Kaur, A. Sharma, and A. K. Sharma, Impact of fear effect on plankton-fish system dynamics incorporating zooplankton refuge, Chaos, Solitons & Fractals 143 (2021), 110563. 12. L. Lai, Z. Zhu, and F. Chen, Stability and bifurcation in a predator–prey model with the additive allee effect and the fear effect, Mathematics 8 (2020), no. 8, 1280. 13. J. Liu, Qualitative analysis of a diffusive predator–prey model with allee and fear effects, International Journal of Biomathematics 14 (2021), no. 06, 2150037. 14. J. Ma, Y. Xu, C. Wang, and W. Jin, Pattern selection and self-organization induced by random boundary initial values in a neuronal network, Physica A: Statistical Mechanics and Its Applications 461 (2016), 586–594. 15. K. Manna and M. Banerjee, Stationary, non-stationary and invasive patterns for a prey-predator system with additive allee effect in prey growth, Ecological complexity 36 (2018), 206–217. 16. K. Sarkar and S. Khajanchi, Impact of fear effect on the growth of prey in a predator-prey interaction model, Ecological Complexity 42 (2020), 100826. 17. S. K. Sasmal and Y. Takeuchi, Dynamics of a predator-prey system with fear and group defense, Journal of Mathe- matical Analysis and Applications 481 (2020), no. 1, 123471. 18. D. Sen, S. Ghorai, S. Sharma, and M. Banerjee, Allee effect in prey’s growth reduces the dynamical complexity in prey-predator model with generalist predator, Applied Mathematical Modelling 91 (2021), 768–790. 19. M. Sen and M. Banerjee, Rich global dynamics in a prey–predator model with allee effect and density dependent death rate of predator, International Journal of Bifurcation and Chaos 25 (2015), no. 03, 1530007. 20. R. Shi and Z. Hu, Dynamics of three-species food chain model with two delays of fear, Chinese Journal of Physics 77 (2022), 678–698. 21. G. Sun, Z. Jin, L. Li, and B. Li, Self-organized wave pattern in a predator-prey model, Nonlinear Dynamics 60 (2010), no. 3, 265–275. 22. G. Sun, L. Li, Z. Jin, Z. Zhang, and T. Zhou, Pattern dynamics in a spatial predator-prey system with allee effect, Abstract and applied analysis, vol. 2013, 2013. 23. X. Tian and S. Guo, Spatio-temporal patterns of predator–prey model with allee effect and constant stocking rate for predator, International Journal of Bifurcation and Chaos 31 (2021), no. 16, 2150249. 24. A. M. Turing, The chemical basis of morphogenesis, Bulletin of mathematical biology 52 (1990), no. 1, 153–197. 25. K. Vishwakarma and M. Sen, Influence of allee effect in prey and hunting cooperation in predator with holling type-iii functional response, Journal of Applied Mathematics and Computing (2021), 1–21. 26. C. Wang, M. Lv, A. Alsaedi, and J. Ma, Synchronization stability and pattern selection in a memristive neuronal network, Chaos: An Interdisciplinary Journal of Nonlinear Science 27 (2017), no. 11, 113108. 27. C. Wang and J. Ma, A review and guidance for pattern selection in spatiotemporal system, International Journal of Modern Physics B 32 (2018), no. 06, 1830003. 28. X. Wang, Y. Tan, Y. Cai, and W. Wang, Impact of the fear effect on the stability and bifurcation of a leslie–gower predator–prey model, International Journal of Bifurcation and Chaos 30 (2020), no. 14, 2050210. 29. L. Wolpert, Positional information and the spatial pattern of cellular differentiation, Journal of theoretical biology 25 (1969), no. 1, 1–47. 30. L. Wolpert, Development of pattern and form in animals, Carolina Biological Supply Co., 1977. 31. Z. Xiao, Z. Li, et al., Stability analysis of a mutual interference predator-prey model with the fear effect, Journal of Applied Science and Engineering 22 (2019), no. 2, 205–211. 32. B. Xie, Impact of the fear and allee effect on a holling type ii prey–predator model, Advances in Difference Equations 2021 (2021), no. 1, 1–15. 33. B. Xie and N. Zhang, Influence of fear effect on a holling type iii prey-predator system with the prey refuge, AIMS Mathematics 7 (2022), no. 2, 1811–1830. 34. X. Yang, H. Qiu, and T. Zhou, Oscillatory turing pattern formation from the interactions between hopf and turing bifurcations, Chinese Journal of Physics (2015). 35. P. Ye and D. Wu, Impacts of strong allee effect and hunting cooperation for a leslie-gower predator-prey system, Chinese Journal of Physics 68 (2020), 49–64. 36. L. Y. Zanette, A. F. White, M. C. Allen, and M. Clinchy, Perceived predation risk reduces the number of offspring songbirds produce per year, Science 334 (2011), no. 6061, 1398–1401. DIFFUSION-DRIVEN INSTABILITY AND PATTERN FORMATION IN A PREY-PREDATOR SYSTEM 99 37. H. Zhang, Y. Cai, S. Fu, and W. Wang, Impact of the fear effect in a prey-predator model incorporating a prey refuge, Applied Mathematics and Computation 356 (2019), 328–337. 38. X. Zhang, Q. An, and L. Wang, Spatiotemporal dynamics of a delayed diffusive ratio-dependent predator–prey model with fear effect, Nonlinear Dynamics 105 (2021), no. 4, 3775–3790. 39. X. Zhang, H. Zhao, and Y. Yuan, Impact of discontinuous harvesting on a diffusive predator–prey model with fear and allee effect, Zeitschrift für angewandte Mathematik und Physik 73 (2022), no. 4, 1–29. 40. Z. Zhu, R. Wu, L. Lai, and X. Yu, The influence of fear effect to the lotka–volterra predator–prey system with predator has other food resource, Advances in Difference Equations 2020 (2020), no. 1, 1–13. D. Pal, Centre for Mathematical Biology and Ecology, Department of Mathematics, Jadavpur University, Kolkata-700032, India Email address: debjitpal381@gmail.com D. Kesh, corresponding author, Centre for Mathematical Biology and Ecology, Department of Mathe- matics, Jadavpur University, Kolkata-700032, India Email address: dkeshju@gmail.com D. Mukherjee, Department of Mathematics, Vivekananda College, Thakurpukur, Kolkata-700063, India Email address: mukherjee1961@gmail.com 1. Introduction 2. Equilibrium Points 3. Local Stability Analysis of Equilibrium Points: 4. Hopf Bifurcation 5. Numerical Simulation 6. Model With Diffusion 6.1. Diffusion-driven Instability 6.2. Numerical simulation 7. Discussion References