ANALYTICAL DETECTION OF STATIONARY TURING PATTERN IN A PREDATOR-PREY SYSTEM WITH GENERALIST PREDATOR

. A prey-predator model with Holling type-II functional response and a generalist predator exhibits complex dynamics in response to parameter variation. Generalist predators implicitly exploit-ing multiple food resources reduce predation pressure on their focal prey species that causes it to become more stable compared to a prey-predator system with specialist predator. In the temporal system, bistability and tristability are observed along with various global and local bifurcations. Existence of homogeneous and heterogeneous positive steady state solutions are shown to exist for suitable ranges of parameter values in the corresponding spatio-temporal diﬀusive system. Weakly nonlinear analysis, using multi-scale perturbation technique, is employed to derive amplitude equation for the stationary patterns near the Turing bifurcation threshold. The analytical results of the amplitude equations are validated using exhaustive numerical simulations. We also identify bifurcation of multiple stable stationary patch solutions as well as dynamic pattern solution for parameter values in the Turing and Turing-Hopf regions.

first showed the significance of spatial heterogeneity in long term survival of species paramecium and didinium through a laboratory experiment. Experimenting on two mite species, Typhlodromus occidentalis and Eotetranychus sexmaculatus, Huffaker [18] observed that both species can survive for long time in appropriate heterogeneous environment whereas homogeneous environment set up can lead to extinction of both species.
Turing and Turing-Hopf bifurcations play significant roles in the emergence of various spatio-temporal patterns [5,12,50]. Spatio-temporal patterns can be classified in two categories: stationary patterns and dynamic patterns. Turing bifurcation is responsible for the formation of stationary heterogeneous patterns through Turing instability which occurs when a temporally stable homogeneous solution becomes unstable due to small amplitude of spatially heterogeneous perturbation and the resulting patterns are known as Turing patterns [25,40,45]. These patterns may resemble like spots (hot and cold), stripes or a mixture of them [2,9,17]. Most of the previous works mainly focused on the dependence of Turing pattern formation on control parameter(s) of the RD system [2,5,50]. However, very little attention has been paid on the mode selection problem of the Turing patterns [1,7,29]. Dynamic spatial patterns, which are usually observed in Turing-Hopf domain [6,39,50], are characterized by periodic, quasi-periodic or even spatio-temporal chaotic patterns with the advancement of time. Apart from these, dynamic spatially heterogeneous patterns are also observed in RD systems due to other mechanisms that include travelling wave, periodic travelling wave, modulated travelling wave and invasion [8,[41][42][43].
Predator-prey models for two species interaction are building blocks for several food-chains and food-webs. In such models, predators can be categorized into two groups, namely, specialist predators and generalist predators. Specialist predators become extinct in absence of prey since it completely depends on prey species for food. In the prey-predator system with specialist predator, persistence of both species is observed in a small range of parameter values and the system is driven to extinction outside the parameter range due to the disappearance of co-existing steady state from the interior of first quadrant. Generalist predators are diverse and abundant and they have significant role in natural systems. Generalist predators can survive in absence of their focal prey as they have alternative food sources other than its focal prey. Compared to the specialist predators, generalist predators have received less attention. Spencer and Collie [44] added a linear growth term with the quadratic death rate in generalist predator population so that they have logistic type growth for the predator in the absence of prey population. Alternative food sources of the generalist predator may include a variety of prey sources which can be modelled by two or more prey population along with the predator population [23]. Classical RD prey-predator models with specialist predator having linear death rate and prey's density dependent numerical response are unable to produce stationary patterns [3,39]. But consideration of intraspecies competition among predators in the same model supports the formation of Turing patterns [4,31]. In addition, prey-predator models with linear death rate of specialist predator and predator's density dependent functional response can support Turing pattern formation too [49,50]. Liao et al. [26] studied the effect of diffusion in a ratio-dependent prey-predator model with generalist predator.
In this work, we have chosen a generalist predator with logistic growth rate and prey-predator interaction is modelled with prey dependent functional response. We first carry out a complete bifurcation analysis of the temporal model and illustrate some representative dynamics for certain range of parameter values. Existence of homogeneous and heterogeneous positive steady state solutions are established for suitable ranges of parameter values in the RD system. Analytical expressions of one dimensional Turing patterns corresponding to unstable eigen modes near the mode selecting Turing thresholds are derived using weakly nonlinear analysis and these are validated with the help of numerical simulations. This paper is organized as follows. In Section 2, we describe our temporal system, existence of its equilibrium points and their stability, followed by Section 3 that analyzes all possible local bifurcations and a parametric bifurcation diagram in Section 4. Next, we consider the corresponding one dimensional spatial extension and prove the existence of global solution and stationary coexistence solution in Section 5. We also check the stability of the homogeneous steady state, derive the Turing instability criteria and various mode selecting thresholds for parameter values in the Turing domain. Amplitude equation is derived using weakly nonlinear analysis in Section 6 and the results of the weakly nonlinear analysis are compared with the numerical results in Section 7. Finally, discussions and conclusions are drawn in Section 8.

Temporal model
We consider a prey-predator model with generalist predator where both the species have intra-specific competition and their interaction is modelled by Holling type-II functional response [36]. Let N (T ) and P (T ) be the prey and predator population density at time T respectively. The system is governed by two ordinary differential equations subject to initial conditions N (0), P (0) ≥ 0. Parameters r 1 , r 2 represent the intrinsic growth rates of the prey and predator population respectively and b 1 , b 2 are the corresponding intra-species competition. The parameter m represents the consumption rate of prey by the predator and h is the handling time, i.e., the average time spent by predator on a prey. Also, e is the conversion coefficient which represents the fraction of consumed prey helping to the predator growth. All the parameters in the system (2.1) are assumed to be positive. We take t = r 1 T, u = b 1 N /r 1 and v = b 2 P /r 2 for dimensionless time, prey population density and predator population density respectively. Using these variables in system (2.1), we obtain and Here, a = b 1 r 2 /b 2 r 2 1 h, α = b 1 /mhr 1 , η = r 2 /r 1 and b = eb 1 r 2 /b 2 r 2 1 h are positive dimensionless parameters.

Positivity and boundedness
Any solution (u(t), v(t)) of system (2.2), starting from a non-negative initial condition (u(0), v(0)), remains non-negative and bounded for all time t. Proof of positivity and boundedness are presented in Appendix A.

Equilibrium points and their existence
The system (2.2) always has trivial equilibrium E 0 (0, 0), and two axial equilibria E 1 (1, 0) and E 2 (0, 1). Intersection of non-trivial prey nullcline f 1 (u, v) = 0 and predator nullcline f 2 (u, v) = 0 in the first quadrant generates co-existing equilibrium points E * j (u j , v j ), where j can vary from 1 to 3. Figure 1. Location of the co-existing equilibrium points for different configurations of the nontrivial nullclines depending on parametric conditions. Here red and green color curves represent the nontrivial predator and prey nullclines respectively.
In the first subcase, the system has three co-existing equilibrium points [see Fig. 1b] and only one equilibrium point in the second subcase [see Fig. 1c] and third subcase [see Fig. 1d]. III. For α < a, we have two sub-cases: (a) η < η SN I ; or (b) η > η SN I . In the first subcase, no interior equilibrium point exists [see Fig. 1e] and two interior equilibria exist in the second subcase [see Fig. 1f].
The values η SN I , η SN II of η and α CP 1 of α are various bifurcation threshold values that we will describe later.

Linear stability of equilibria
To find the linear stability of an equilibrium point E * (u * , v * ), we linearize the system (2.2) about the equilibrium point. This leads to the Jacobian matrix of the system (2.2) evaluated at E * (u * , v * ) given by The equilibrium point E * is stable (resp. unstable) if both the eigenvalues of J(E * ) have negative (resp. positive) real part. On the other hand, the equilibrium point is called a saddle point when the eigenvalues have opposite sign.
Proposition 2.1. The following statements hold for the system (2.2): (i) The trivial fixed point E 0 is always unstable.
(ii) The axial fixed point E 1 is always a saddle point. But the other axial fixed point E 2 is asymptotically stable for a > α and a saddle fixed point for a < α.
Proof. (i) The eigenvalues of J(E 0 ) are 1 and η and hence E 0 is always unstable.
(ii) The Jacobian matrix calculated at E 1 has one positive eigenvalue η + b 1+α and one negative eigenvalue −1. Hence, E 1 is a saddle point. The Jacobian calculated at E 2 is given by Thus, E 2 is a saddle point for α > a with stable manifold along the v axis. On the other hand, E 2 becomes asymptotically stable for α < a.
It is difficult to determine the nature of co-existing equilibrium points such as E * 1 , E * 2 and E * 3 (whenever they exist) analytically due to the unavailability of explicit expression of these points. Existence of these equilibria and their stability in different regions depending on parameter values will be discussed in Section 4.
Proof. From (2.7), the Jacobian matrix J(E 2 ) at α = α T C has a zero eigenvalue. Let the eigenvectors corresponding to the zero eigenvalue of the Jacobian matrix and its transpose be p = [ ηa b , 1] T and q = [1, 0] T respectively. Then transversality conditions [38] become Here, F = [F 1 (u, v), F 2 (u, v)] T and all other notations are same as in [38]. Thus, the system (2.2) satisfies all the transversality conditions of transcritical bifurcation at E 2 for α = α T C .

Saddle-node bifurcation
Whenever a 2 < 0, a 1 > 0 and a 0 < 0 hold, then the system (2.2) has either three co-existing equilibria or a unique co-existing equilibrium point. Suppose u s1 and u s2 are two real roots of G (u) = 0 with 0 < u s1 < u s2 . If G(u s1 ) and G(u s2 ) are of opposite sign, then three co-existing equilibria exist and we denote them by E * 1 , E * 2 and E * 3 with u 1 < u 2 < u 3 . Now two cases may arise due to variation of one temporal parameter: either E * 2 coincides with E * 1 when G (u s1 ) = 0 or E * 2 coincides with E * 3 when G (u s2 ) = 0. Thus, two saddle-node bifurcations SN I and SN II occur for the first case and second case, respectively. We choose η as a saddle-node bifurcation parameter and the threshold value for SN I is given by For the SN II threshold value, u s1 is replaced by u s2 in the right side of the expression for η SN I . The transversality conditions for this bifurcation is shown in the following proposition.
Therefore, one of the eigenvalues of the Jacobian matrix J(E * SN ) must be zero. Let p = 1, (3.1b) Hence, the system (2.2) undergoes saddle-node bifurcation at η = η SN .

Cusp bifurcation
We have observed that two saddle-node bifurcation SN I and SN II occur at E * SN I = (u s1 , v s1 ) and E * SN II = (u s2 , v s2 ) for varying the bifurcation parameter η. If we alter another temporal parameter α, then E * SN I and E * SN II coincide at E * CP1 = (u CP1 , v CP1 ) and these two saddle-node bifurcation curves SN I and SN II intersect at a cusp bifurcation point (η CP 1 , α CP 1 ) in the η − α parametric domain. Proof. Suppose at η = η CP 1 and α = α CP 1 , G(u) has triple root u CP1 , i.e., G(u CP1 ) = G (u CP1 ) = G (u CP1 ) = 0, where Feasibility condition for E * CP1 = (u CP1 , v CP1 ) comes from a 2 < 0. By following the same procedure as in Proposition 3.2, one can show that the Jacobian evaluated at E * CP1 has a zero eigenvalue. Let an eigenvector of J(E * CP1 ) and [J(E * CP1 )] T corresponding to zero eigenvalue be p and q respectively. Then, the transversality condition (3.1b) for the cusp bifurcation becomes which shows that the system (2.2) undergoes a codimension-2 cusp bifurcation at the threshold value (η CP 1 , α CP 1 ). Two saddle-node bifurcation curves SN I and SN II emerge from this cusp bifurcation point. The relationship between η CP 1 and α CP 1 is given by Proposition 3.4. The system (2.2) undergoes another cusp bifurcation at E 2 when G(0) = G (0) = 0.

Hopf bifurcation
A stable interior equilibrium E * can lose its stability via Hopf Bifurcation when the trace of the Jacobian matrix evaluated at E * changes from negative to positive value due to variation of a model parameter. Here, we choose α as a Hopf bifurcation parameter. The Jacobian matrix calculated at E * (u * , v * ) is given by Now let us assume that α = 2) undergoes a Hopf bifurcation at α = α H if the following non-hyperbolicity and transversality conditions are satisfied: Proposition 3.5. Interior equilibrium point E * 1 of the system (2.2) becomes unstable via Hopf bifurcation at α = α H .
A Hopf-bifurcation is called supercritical (resp. subcritical) if the limit cycle generated at the Hopf-bifurcation is stable (resp. unstable). Explicit expression of Hopf bifurcation threshold is hard to find due to the lack of explicit expression of the co-existing fixed point. But we can verify these results numerically. If we fix the temporal parameters a = 0.2, b = 0.15 and η = 0.095, then the system (2.2) undergoes subcritical Hopf bifurcation around E * 1 = (0.04801, 1.28156) for α H = 0.22122. For the same set of parameter values except η = 0.0915, supercritical Hopf bifurcation observed around E * 1 = (0.05084, 1.30333) for α H = 0.22379.

Bautin bifurcation
We have seen that the co-existing equilibrium point E * 1 changes its stability either through a supercritical or a subcritical Hopf bifurcation. A supercritical Hopf corresponds to the first Lyapunov coefficient l 1 < 0 in which a stable limit cycle is formed. On the other hand, first Lyapunov coefficient l 1 is positive in subcritical Hopf bifurcation and an unstable limit cycle arises. In case of zero first Lyapunov coefficient, the system (2.2) passes through a codimension-2 bifurcation called a generalized Hopf bifurcation (GH) or a Bautin bifurcation.
Analytical expression and transversality conditions for this bifurcation are hard to find due to the same reason mentioned before. But numerically we verify that for parameter values a = 0.2 and b = 0.15, the system (2.2) undergoes two Bautin bifurcation at α GH 1 = 0.23412 and η GH 1 = 0.09399 and at α GH 2 = 0.22049 and η GH 2 = 0.09705.

Bogdanov-Takens bifurcation
We have observed that system (2.2) has two saddle-node bifurcation curves and two Hopf bifurcation curves in the η − α parametric plane. When one Hopf bifurcation curve intersects with a saddle-node bifurcation curve then a codimension-2 bifurcation occurs. This bifurcation is called Bogdanov-Takens (BT) bifurcation. Generally, BT bifurcation occurs at a fixed point when both the trace and the determinant of the Jacobian matrix evaluated at that point vanish.
Proposition 3.7. The system (2.2) undergoes a Bogdanov-Takens Bifurcation at E * 1 or E * 3 if the Jacobian matrix at E * 1 or E * 3 has zero eigenvalue of multiplicity two. Again, it is difficult to find analytical expressions and corresponding transversality conditions for the BT bifurcation since explicit expressions for co-existing equilibria are not available. However, we choose η and α as bifurcation parameters and verify the existence of BT bifurcation numerically for the system (2.2). For this, we fix other parameters at a = 0.2, b = 0.15 and find that the system (2.2) undergoes a BT bifurcation for E * 1 at (η BT 1 , α BT 1 ) = (0.10037, 0.22556) and another BT bifurcation for E * 3 at (η BT 2 , α BT 2 ) = (0.12413, 0.21001).

Bifurcation diagram
Here we discuss the global and local bifurcation results of the system (2.2) obtained with the help of Matlab software package. We present a bifurcation diagram in the two-dimensional η − α plane for fixed values of a = 0.2 and b = 0.15. The bifurcation threshold values are very close to each other. For example, if we take η = 0.1, then the values of other bifurcation thresholds are α SN I = 0.22604, α SN II = 0.226098, α H 1 = 0.226044 and α SN LC I = 0.226184. These values won't be clearly visible in actual bifurcation plot. Hence, a schematic bifurcation diagram is plotted in Figure 2 for better visualisation. Here, α SN LC I denotes the threshold value of saddle node bifurcation of limit cycle at which an unstable limit cycle and a stable limit cycle surrounding an interior equilibrium point coincide. Figure 2. Schematic bifurcation diagram in η − α plane. Cyan, blue and red color curves respectively represent transcritical, saddle-node and Hopf bifurcation curves. Green color curve from GH point represents saddle-node bifurcation of limit cycle and magenta color curve from BT point represents homoclinic orbit. Another global homoclinic bifurcation curve is represented by brown color curve.
There are local as well as global bifurcation curves. Local bifurcation curves include saddle-node bifurcation curve (blue color curve), transcritical bifurcation curve (cyan color curve) and Hopf bifurcation curve (red color curve). On the other hand global bifurcation curves include saddle-node bifurcation of limit cycle (green color curve) and homoclinic curve (magenta and brown color curves). These curves divide the η − α plane into fourteen regions R 1 − R 14 . From Proposition 3.2, the trivial equilibrium E 0 is unstable and axial equilibrium point E 1 is a saddle for all the regions. However, the other axial equilibrium point E 2 is stable in regions R 5 and R 6 , and a saddle point elsewhere. Now we discuss the behavior of the interior equilibrium points of the system (2.2) in each region and the changes in dynamical behavior across the bifurcation curves. Dynamics of each region is listed in Table 1 and schematic phase portraits (for clarity) of the co-existing equilibria in each region are plotted in Figure 3. It is seen that the regions R 1 , R 7 have similar dynamics and so are the regions R 11 , R 13 .
We first discuss local bifurcation curves across which number of equilibria or stability change with the variation of parameters. Figure 2 shows that two saddle-node bifurcation curves emerge from a cusp bifurcation point CP 1 (η CP 1 , α CP 1 ) ≡ (0.099, 0.227) which is a codimension-2 bifurcation point. We denote the lower and upper saddle-node bifurcation curves by SN I and SN II respectively. The system (2.2) has three interior equilibria E * 1 , E * 2 and E * 3 in the region bounded by SN I , SN II and the transcritical bifurcation curve (T C) parallel to η-axis. We label the three equilibria E * 1 , E * 2 and E * 3 such that their u-components satisfy 0 < u 1 < u 2 < u 3 . Among them E * 2 is always saddle and the nature of E * 1 and E * 3 depend on the Hopf bifurcation curve. In region R 2 , E * 1 is an unstable node, E * 2 is saddle and E * 3 is a stable node. E * 1 and E * 2 collide and disappear across the SN II curve leading to a unique equilibrium point. Whenever the system (2.2) has only a unique co-existing equilibrium point then we label it as E * 1 . Thus, the system has unique co-existing equilibrium point E * 1 in R 1 , which is globally asymptotically stable. Hopf curves meet with saddle-node bifurcation curves at codimension-2 BT bifurcation points BT 1 and BT 2 . We denote Hopf curves emanating from BT 1 and BT 2 by H 1 and H 2 respectively. Stability of co-existing equilibrium point changes when α crosses the Hopf threshold value α H 1 or α H 2 . There are another two codimension-2 Bautin (or GH) bifurcation points GH 1 and GH 2 on the Hopf curves H 1 and H 2 respectively, where the first Lyapunov co-efficient (l 1 ) vanishes. A sub-critical (resp. super-critical) Hopf bifurcation occurs across each Hopf curve to the right (resp. left) of the GH point. Numerical threshold Table 1. Dynamical behaviour of the interior equilibria in each region for the system (2.2). Axial fixed point E 2 is saddle in all regions except in R 5 and R 6 where it is stable.

Interior equilibria Remarks
There is a bistablity between E * 1 and E * 3 .
There is a bistablity between axial fixed point E 2 and interior fixed point E * 3 . R 6 (Fig. 3f) ---Axial equilibrium E 2 is asymptotically stable. R 7 (Fig. 3a) Stable ---R 8 (Fig. 3g) Unstable --A stable Hopf bifurcating limit cycle appears around the co-existing fixed point E * 1 . R 9 (Fig. 3h) Unstable Saddle Unstable All the interior equilibrium points lie inside a stable limit cycle.
3 are surrounded by a stable limit cycle and the system shows bistable phenomena where the stable limit cycle and a stable co-existing equilibrium point act as attractor. R 11 (Fig. 3j) Stable --E * 1 is encircled by an unstable limit cycle which is again encircled by a stable limit cycle. R 12 (Fig. 3k) Stable Saddle Stable An unstable limit cycle containing E * 1 , E * 2 and E * 3 are surrounded by a stable limit cycle. R 13 (Fig. 3j) Stable --Same characteristic as in R 11 . R 14 (Fig. 3l) Unstable Saddle Stable A stable limit cycle surrounds all the co-existing fixed points.
values of the GH and BT bifurcation points have already been given. T C curve intersects with SN II curve at another cusp bifurcation point CP 2 (0.2,0.1875). Interior equilibrium E * 1 and axial equilibrium E 2 exchange their stability across T C. However, E * 1 becomes infeasible in region R 5 and hence the system (2.2) has two co-existing equilibria E * 2 (saddle) and E * 3 (stable node) in R 5 . The stable manifolds of E * 2 separates the basins of attraction of E 2 and E * 3 . E * 2 and E * 3 collide through saddle-node bifurcation (across SN I ) and disappear in region R 6 . Thus, the system (2.2) has no co-existing equilibrium point and E 2 is globally asymptotically stable in R 6 . From R 6 to R 7 across T C curve, axial point E 2 becomes unstable and E * 1 again becomes feasible. Hence, E * 1 is globally asymptotically stable in region R 7 . Now we discuss global bifurcations which occur when there is a collision between two periodic orbits or a periodic orbit and an equilibrium point. A saddle-node bifurcation curve of limit cycles (SN LC I ) emerges from the GH 1 point and meets with the SN II curve. On the SN LC I curve, an unstable limit cycle surrounding E * 1 generated from subcritical Hopf bifurcation coincides with a stable limit cycle. Therefore, stable equilibrium point E * 1 is surrounded by an unstable limit cycle inside a stable limit cycle in R 13 . A homoclinic curve (HC I ) emerges from the BT 1 point and ends at the intersection point of SN LC I and SN II curves. HC 1 corresponds to the appearance or disappearance of an unstable limit cycle around the co-existing fixed point E * 1 . A second saddle-node bifurcation curve of limit cycles (SN LC II ) from the GH 2 point and a homoclinic curve (HC II ) from the BT 2 point meet at a point on the SN I curve.
Another homoclinic curve HC III (brown color curve) extends from the intersection point of HC I , SN LC I and SN II to the intersection point of HC II , SN LC II and SN I . Along the HC III curve, a large stable limit cycle containing all the co-existing equilibria E * 1 , E * 2 and E * 3 appears or disappears. There is bistablity between Figure 3. Schematic phase portrait diagram. Here green circle represents stable limit cycle and red circle represents unstable limit cycle. Stable and unstable manifolds of saddle point are denoted by cyan and magenta color curves respectively.
the stable limit cycle and the co-existing fixed point E * 3 in R 10 and the stable manifold of the saddle point E * 2 separates their basin of attraction.

Spatio-temporal model
Now we take the spatial effects into account by including the diffusion terms. For simplicity, we consider the spatio-temporal model in one dimension. Let Ω = (0, L) be a bounded domain with boundary ∂Ω = {0, L}. The governing equations are where d represents the ratio of diffusion coefficients of predator and prey. Next we discuss existence of global solution in Section 5.1, existence of positive steady state solution in Section 5.2 and stability of homogeneous steady states in Section 5.3. The results obtained are applicable in higher dimensions as well. Hence, we take ∇ 2 for Laplacian and ∂ ∂n for normal derivative in Sections 5.1-5.3.

Existence of global solution
The local existence of unique solution of the system (5.1) is easy to establish by standard existence theory [14,22]. Here, we prove the global existence by showing the boundedness of ||u(t, ·)|| L∞ , ||v(t, ·)|| L∞ for all finite time t. Let D T = (0, T ] × Ω and S T = (0, T ] × ∂Ω. To prove the global existence, we make use of the following strong comparison theorem [43,48]: t, x), . . . , p m (t, x)) T and q(t, x) = (q 1 (t, x), . . . , q m (t, x)) T , m ≥ 1, be two vector functions that satisfy the following conditions (where inequalities satisfy component-wise): Here, each f i is Hölder continuous for (t, x) ∈D T and continuously differentiable with respect to each of the components p 1 , p 2 , . . . , p m . Further, each f i is quasi-monotonically increasing in p, i.e., Then, p(t, x) ≤ q(t, x) inD T .
Proposition 5.2. For any non-negative initial condition u 0 (x), v 0 (x) ∈ C 2 (Ω) with the system (5.1) has unique global and bounded solution for all t ≥ 0.
Proof. Using Lemma 14.20 of [43], it can be shown that u(t, x) ≥ 0 and v(t, x) ≥ 0 for any non-negative initial condition u 0 (x), v 0 (x). Let A(t) = A0 A0+(1−A0)e −t be the solution of the initial value problem Consider the following problem: Then, using Theorem 5.1, we conclude that u(t, x) ≤ A 0 ≤ 1 and hence u(t, x) ≤ 1 inD T . To show the boundedness of v, we consider the initial value problem The solution of the system (5.2) is From second equation of system (5.1), we obtain Then, using Theorem 5.1, we find v(t, x) ≤ B 0 and hence v(t, x) ≤ (η+b) η inD T . Thus, ||u(t, .)|| L∞ , ||v(t, .)|| L∞ are bounded for all finite t.

Existence of positive steady state solutions
A steady state solution (u(x), v(x)) of the system (5.1) satisfy  Proof. For u ≥ 0 and v ≥ 0, Hence, the system (5.3) is a mixed quasi-monotone [37,43]. Now, an upper solution (u(x), v(x)) and lower solution (u(x), v(x)) of the system (5.3) satisfy the following: To prove the existence of positive steady state, we construct a pair of upper and lower solutions (u(x), v(x)) and (u(x), v(x)) of the system ( We also require v(x) = M to satisfy the right hand side of the inequality (5.4), i.e., For a < 1 and α − a + 1 − M > 0, it is easy to verify that u(x) = 1 − a satisfies the inequality (5.6). For the given u(x) = 1 − a with a < 1 and 1 − a + α > M , if we choose v(x) to be a small positive constant such that then the right hand side of the inequality (5.5) is satisfied. Thus, we have constructed a pair of positive upper and lower solutions (u(x), v(x)) and (u(x), v(x)) of the system (5.3). This completes the proof.
Proof. Linearized system (5.7) around E 0 becomes Now E 0 is unstable if the largest eigenvalue of the corresponding eigenvalue problem is positive. Let λ m be the largest eigenvalue of the system (5.8). Note that the principal eigenvalue λ 1 of the system is positive with corresponding eigenfunctionw 1 . But, λ 1 is also an eigenvalue of the system (5.8) with the eigenfunction (w 1 , 0) T . Therefore, λ m ≥ λ 1 > 0. This completes the proof.
E 1 is unstable if the largest eigenvalue λ = λ m of the corresponding eigenvalue problem is positive. Now the principal eigenvalue λ 1 of is positive with corresponding eigenfunctionw 2 . Letw 1 be the solution of the linear problem Then, λ 1 is also an eigenvalue of the system (5.9) with the eigenfunction (w 1 ,w 2 ) T . Therefore, λ m ≥ λ 1 > 0. Hence the proof is complete.
Corresponding eigenvalue problem is (5.10) If the largest eigenvalue λ m of the system (5.10) is positive then E 2 is unstable and stable otherwise. The principal eigenvalue λ 1 of the system is positive for α > a with corresponding eigenfunctionw 1 . Suppose thatw 2 is the solution of the linear problem Then, λ 1 is also an eigenvalue of the system (5.10) with the eigenfunction (w 1 ,w 2 ) T . Therefore, λ m ≥ λ 1 > 0.
Hence, E 2 is unstable for α > a. For α < a, let (w 1 ,w 2 ) T be the eigenfunction of the system (5.10) corresponding to the largest eigenvalue λ m . Ifw 1 = 0, then λ m is also an eigenvalue of the system (5.11) corresponding to eigenfunctionw 1 . Now α < a implies λ m < 0. Ifw 1 = 0, thenw 2 = 0 and λ m becomes an eigenvalue with eigenfunctionw 2 for the problem Clearly, the largest eigenvalue of the system (5.12) is −η < 0 and hence λ m < 0. As a result, the steady state solution E 2 (0, 1) is stable for α < a.
Proof. The linearized system about E * is ∂W ∂t = D∇ 2 W + JW, (5.13) where D =diag(1, d), W is a two-dimensional column vector and Let 0 = k 1 < k 2 < · · · < k j < · · · be the eigenvalues and E(k j ) be the eigenfunction space corresponding to k j for the eigenvalue problem Further, suppose that {ψ i,j : i = 1, . . . , dim E(k j ) } be the orthogonal basis set of E(k j ) and We expand the solution W as where C j (t) is a two-dimensional column vector and ψ j (x) ∈ Ψ. Substituting (5.14) into (5.13) and equating the coefficient of each ψ j , we find is also positive for a 10 < 0 since a 01 < 0 and b 10 > 0. Now a 10 < 0 implies (u * +α) 2 v * > a. Thus, steady state

Turing instability
We consider the parameter restriction η < η CP 1 and α > a so that the local model (5.1) has unique co-existing homogeneous steady state E * (u * , v * ). For Turing instabilities, the homogeneous steady state solution is stable in the absence of diffusion, i.e., a 10 + b 01 < 0 and (a 10 b 01 − a 01 b 10 ) > 0. To check Turing instability, we introduce spatial perturbation around E * by u = u * + ūe (λt+ikx) and v = v * + ve (λt+ikx) , where | | << 1. Substituting u and v into the system (5.1) and linearizing we find For non-trivial solution of the system (5.15) we have det(S k ) = 0, which leads to For instability, we must have either T (k 2 ) > 0 or H(k 2 ) < 0 hold for some k. We observe that T (k 2 ) < 0 holds for all k since a 10 + b 01 < 0. Hence, the condition for Turing instability is H(k 2 ) < 0 for some k. The minimum value of H(k 2 ) is Hence, we must have H min = 0 at the Turing bifurcation threshold and the corresponding critical value of d is given by .
For d > d c , the system (5.1) shows non-homogeneous stationary pattern of certain mode under spatially heterogeneous perturbation.

Mode selection
Suppose that for d > d c , there exists some k ∈ (k L , k R ) such that H(k 2 ) < 0, where Note that k 2 L and k 2 R are obtained from solving H(k 2 ) = 0. If we take Ω = [0, L], then the n-th mode wave number is k n = nπ/L, where n is a natural number. If k n lies between k L and k R , then it is an unstable mode and the system (5.1) can show n-th mode pattern solution. This n-th mode selecting bifurcation threshold value is given by Considering η as a bifurcation parameter, we find the corresponding value of d (n) c and denote this n-th mode selecting threshold curve [34,35] by M n in the η − d plane. The region above the M n curve is unstable with respect to the n-th mode solution (see Fig. 4). Now we take d > d small and investigate the appearance of n-th mode pattern solution of the system (5.1). We also examine the amplitude of the resulting pattern solution.

Weakly nonlinear analysis
In this section, we apply weakly nonlinear analysis using multi-scale perturbation method and Fredholm theory to derive amplitude equations [7,16,29] for the n-th mode Turing patterns solution arising near the critical c ). First, we expand system (5.1) around E * (u * , v * ) using Taylor series expansion upto fifth order: Dropping the tilde for convenience, we write (6.1) in matrix form as c ), the evolution of pattern amplitude changes slowly. Thus to capture the amplitude correspond to n-th mode wavenumber, we introduce slow temporal multiple scale where is the small perturbation parameter that measures the distance between bifurcation parameter d and the bifurcation threshold (d = d (n) c ), i.e., The solution of the system (6.2) is also expanded in terms of : where P j = (u j , v j ) T . These lead to where s i is the coefficient of i in the expansion of the nonlinear term. Further, Substituting these expression into equation (6.2) and equating like power of , we obtain: where the expressions for F , G, H and I are given in Appendix B. Solution of (6.3) under Neumann boundary condition is and A is the amplitude of the pattern which is still unknown at this stage. Now substitution (6.8) into (6.4) gives rise to a non-homogeneous system. For the solution of such a nonhomogeneous system to exist, the right hand side must be orthogonal to the kernel of the adjoint of L c . This is known as Fredholm solvability condition which imposes d 1 = 0 and T 1 = 0. Substitution of solution of (6.4) into (6.5) and applying solvability condition again, we obtain Stuart-Landau equation for the amplitude A: Equation (6.9) is same as the normal form of pitchfork bifurcation. The detailed derivation of (6.9) is given in Appendix B. Note that σ is always positive but the sign of l depends on G 31 (see Appendix B). We have supercritical case for l > 0 and subcritical case for l < 0.

The supercritical case
The system (6.9) has two stable equilibria A * = ± σ l and an unstable equilibrium A * = 0. These stable equilibria are responsible for pattern formation for d > d (n) c and they attract all non-zero solutions in long time, i.e., Hence the solution of (6.2) is given by c . We validate this result in the next section.

The subcritical Case
When l < 0, then the system (6.9) has unique equilibrium solution A = 0 which cannot represent an amplitude of a pattern. In this case, we need to carry the weakly nonlinear analysis up to O( 5 ). This leads to the following quintic Stuart-Landau equation for the amplitude A (see Appendix C): dA dt =σA −lA 3 +rA 5 , (6.10) whereσ = 2σ + 4 σ ,l = 2l + 4 l andr = 4 r . Ifσ > 0,l < 0 andr < 0 then the equation (6.10) predicts the long-term behavior of the amplitude of the pattern for d > d , i.e., Hence, the solution of (6.2) is given by U = A ∞ Φ cos(k n x) + O( 2 ) at d = d

Numerical Simulations
We fix temporal parameter values a = 0.2, b = 0.15 and α = 0.3 and keep the domain size L = 60 in all the numerical simulations. We have chosen the initial data to be where ε = 10 −2 and n is the desired Turing mode. For η = 0.063 and d = 12, the dispersion relation is shown in Figure 5. We observe that the largest real part of the growth rate λ(k 2 ) is positive for k = k 3 , k 4 and k 5 only. We have also plotted the corresponding M 3 , M 4 and M 5 mode curves in η − d plane together with Turing curve and temporal Hopf curve in Figure 4. If we decrease η from 0.07 to 0.05 with fixed d = 12, then the location of (η, d) moves from stable region to pure Turing region to Turing-Hopf region (see Fig. 4). We plot the corresponding bifurcation diagram (see Fig. 6) in terms of spatial average of u defined by  In this figure, green color curve represents u av for the homogeneous steady state solution which is stable (solid curve) in the stable region and unstable (dotted curve) otherwise. In the Turing region, Turing pattern solution branches of 3-, 4-and 5-modes emerge from η = 0.0643, 0.0647 and 0.0633 respectively. Red color curve denotes the oscillatory pattern solution emerging at the temporal Hopf threshold in the Turing-Hopf domain. There exists bistablity between oscillatory solution and non-homogeneous steady state solutions in the Turing-Hopf domain.
Next we consider d as a bifurcation parameter to study the amplitude of Turing pattern in the pure Turing region for fixed η = 0.064. From (5.17), we see that d The system (7.1) has stable equilibria A ∞ = ±0.1732 √ d 2 and the 4-mode pattern solution is stable. Solution u of the system (5.1) using weakly nonlinear analysis at long time t is given by c ). We compare the numerical solution with that of weakly nonlinear analysis for d = 10.2 in Figure 7 and the two solutions agree well. For 3-mode solution, we find the Stuart-Landau equation of third order as c + 2 d 2 + 4 d 4 , the maximum of the stable stationary pattern solution is given by max(u) = u * + A ∞ . We have plotted the maximum of u from numerical as well as weakly nonlinear analysis in Figure 8. For the numerical continuation of stationary solution brunch, we have used pseudo-arc-length continuation method [19,32]. The system (5.1) shows homogeneous steady state solution represented by blue solid curve for d < d (4) c . Homogeneous steady solution becomes unstable for d > d (4) c , which is marked by blue dotted line. Magenta color curve represents the 4-mode solution branch emerging at d = d (4) c . We denote 3-mode solution brunch by red color curve appearing at d = d (3) c . And black color curves represent both the 3-and 4-mode solutions from weakly nonlinear analysis. Bistablity between the 3-mode solution and 4-mode solution is observed for d > d (3) c . In this case, the system (5.1) does not exhibit hysteresis since the homogeneous steady-state loses its stability through a super-critical steady-state bifurcation.
In order to obtain hysteresis cycle, we choose η = 0.066 for which 3-mode curve is nearest to Turing bifurcation curve. The amplitude equation is given by whereσ > 0,l < 0 andr < 0. The stationary solution of the system (7.4) is plotted in Figure 9a. Homogeneous steady state is stable for d < d c . These unstable branches change direction and become stable at d = d s . Hence, the system (5.1) has stable 3-mode solutions along with homogeneous stable solution for d s < d < d (3) c . Different stable states lead to possible hysteresis as the bifurcation parameter d is varied. Based on periodic variation of the bifurcation parameter d shown in Figure 9a, we demonstrate a hysteresis cycle in Figure 10. Increasing the parameter d above d

Discussion
Modelling prey-predator interaction and understanding their dynamics have important roles in ecology. Most of the works on prey-predator interaction have been carried out with specialist predators. However, many predators have alternative food resources in addition to the focal prey. Hence, we have considered a preypredator model with generalist predator and Holling type-II functional response for predation. The temporal and spatio-temporal dynamics under various parametric conditions are investigated. Various complex dynamics such as bistability, tristability, global bifurcations, etc., and their occurrences with parameter values are explored for the temporal model. Most of the literature on spatio-temporal models deal with Turing threshold and the corresponding weakly nonlinear analysis [1,7,29] around the Turing threshold. In contrast, we have calculated thresholds for various Turing modes and the weakly nonlinear analysis has been carried out around these mode thresholds.
The temporal model exhibits various local bifurcations consisting of saddle-node, transcritical, Hopf, cusp and BT bifurcations. It also shows global bifurcations like homoclinic bifurcation and saddle node bifurcation of limit cycle. Both species persist for α > a either in a co-existing steady state or in an oscillatory state. Note that α > a implies r 1 b 2 > mr 2 . Thus, persistence of both species can occur when intrinsic growth rate of the prey or intra-species competition in predator become large. It also persists when intrinsic growth rate of predator or consumption rate of prey become small. The prey species go to extinction when parameter values lie in region R 6 for any initial population. The generalist predator population survives due to the availability of alternative food resources. This happens for α < a and η < η SN I . Since η = r 2 /r 1 , the conditions for the existence of predator only population is r 1 b 2 < mr 2 and r 2 /r 1 < η SN I . Combining both we find b 2 /m < r 2 /r 1 < η SN I . Thus, the ratio of intrinsic growth rates of the predator and prey must lie between b 2 /m and η SN I for the predator only population to survive. If we increase the value of η beyond η SN I keeping α < a, then we have a bi-stable region (R 5 ) in which both the predator only or a co-existing state persist. In terms of dimensional parameters, this bi-stable region exists for r 2 /r 1 > max{b 2 /m, η SN I }. Hence, the bi-stable region consisting of predator only or a co-existing state exists when the intrinsic growth rate of predator far exceeds the intrinsic growth rate of prey. The system (2.2) admits four different types of bistability: (i) one stable axial equilibrium point and a stable coexistence equilibrium point (region R 5 ), (ii) two stable coexistence equilibria (region R 3 and R 4 ), (iii) a stable coexistence equilibrium point and a stable limit cycle enclosing another unstable coexistence equilibrium point (region R 10 , R 13 and R 14 ) and (iv) a stable equilibrium point surrounded by an unstable limit cycle which is again enclosed by a stable limit cycle (region R 11 ). For the first three cases, the basins of attractions of the attractors are separated by the stable manifold of a saddle coexistence equilibrium point. In the last case, the basins of attraction are separated by an unstable limit cycle. We have also found tri-stability (region R 12 ), where two stable co-existing equilibrium points and a stable limit cycle surrounding them are three attractors. When the system has more than one attractor, the final state of the population species are determined by the initial conditions.
We have also investigated the behaviour of the system in presence of diffusion terms that account for the random movement of the population. We have shown the existence of global solution by showing the boundedness of both the species for all finite time. Also, constructing upper and lower solutions under suitable parametric restriction, we have shown the existence of co-existing non-homogeneous stationary solution. The local stability of the homogeneous steady states have been analyzed with the help of corresponding eigenvalue problems. We have derived the criteria for Turing instability. Further, we have also calculated the threshold values for various Turing mode solutions. The threshold curves representing various unstable Turing modes sometimes cross each other in a two-parametric plane. Also, the unstable Turing mode, which is near the Turing bifurcation threshold, depends on the size of the domain and parameter values. By employing weakly non-linear analysis, we have derived amplitude equation corresponding to n-th unstable mode Turing solution. These solutions corresponding to different modes appear either due to sub-critical or super-critical bifurcation across the mode-thresholds. For sub-critical bifurcation, multiple branches of stable and unstable solutions exist within a range of bifurcation parameter, which may lead to a hysteresis cycle (see Fig. 9a). We have found good agreement between the results of the weakly non-linear analysis and the results of the full numerical solutions of the system (5.1) near the thresholds (see Fig. 8). Multiple stationary Turing patterns corresponding to different unstable modes depend on the choice of initial conditions and these are usually observed in the Turing region. We find the coexistence of stationary and oscillatory solution, depending upon the initial conditions, for parameter values in Turing-hopf domain.
Appendix A. Proof of positivity and boundedness which shows that u(t) ≥ 0, v(t) ≥ 0 for all t ≥ 0 whenever u 0 ≥ 0, v 0 ≥ 0. Hence all solutions starting from non-negative initial conditions remain within the first quadrant of the u-v plane.
To prove the boundedness of the solutions, let us consider two cases: u 0 ≤ 1 and u 0 > 1. For the first case, our claim is that u(t) ≤ 1 for all t ≥ 0. If not, let there be t 1 , t 2 such that u(t 1 ) = 1 and u(t) > 1 for t 1 < t < t 2 . Then from equation (A.1a), we find u(t) = u 0 exp f 1 (u(s), v(s))ds < u(t 1 ) t 1 < t < t 2 , since f 1 (u, v) < 0 for t 1 < t < t 2 . This contradiction implies that u(t) ≤ 1 for all t ≥ 0 whenever u 0 ≤ 1. For the second case, as long as u(t) ≥ 1, we have u(t) = u 0 exp t 0 f 1 (u(s), v(s))ds < u 0 .
Combining both the cases we have u(t) ≤ max{u 0 , 1} for all t ≥ 0. To show the boundedness of v, we consider X(t) = u(t) + v(t). Using (2.2), we find The above inequality leads to X ≤ C − (C − X(0)) exp(−t) and hence X remains bounded as exp(−t) < 1 for t > 0. Hence v(t) is also bounded. Moreover, A = (u, v) ∈ R 2 + : 0 ≤ u ≤ 1, 0 ≤ u + v ≤ C is the positively invariant region for the system (2.2) and all the solutions converge towards the attractor A.