Accepted Manuscript Interaction between prey and mutually interfering predator in prey reserve habitat: pattern formation and the Turing-Hopf bifurcation Jai Prakash Tripathi, Syed Abbas, Gui-Quan Sun, Debaldev Jana, Cui-Hua Wang PII: DOI: Reference: S0016-0032(18)30507-6 https://doi.org/10.1016/j.jfranklin.2018.07.029 FI 3576 To appear in: Journal of the Franklin Institute Received date: Revised date: Accepted date: 4 January 2018 6 June 2018 31 July 2018 Please cite this article as: Jai Prakash Tripathi, Syed Abbas, Gui-Quan Sun, Debaldev Jana, Cui-Hua Wang, Interaction between prey and mutually interfering predator in prey reserve habitat: pattern formation and the Turing-Hopf bifurcation, Journal of the Franklin Institute (2018), doi: https://doi.org/10.1016/j.jfranklin.2018.07.029 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. ACCEPTED MANUSCRIPT Interaction between prey and mutually interfering predator in prey reserve habitat: pattern formation and the Turing–Hopf bifurcation Jai Prakash Tripathi1 , Syed Abbas2 , Gui-Quan Sun3,† , Debaldev Jana4 , Cui-Hua Wang3 1 Department of Mathematics, Central University of Rajasthan, NH-8, Bandarsindri, Kishangarh-305817, Distt.-Ajmer, Rajasthan, India 2 School of Basic Sciences, Indian Institute of Technology Mandi, Mandi-175001, Himachal Pradesh, India Systems Research Center, Shanxi University CR IP T 3 Complex Taiyuan, Shanxi 030006, People’s Republic of China 4 Department of Mathematics & SRM Research Institute, SRM Institute of Science and Technology, Kattankulathur-603 203, Tamil Nadu, India † Corresponding Author. Email: sunguiquan@sxu.edu.cn Abstract ED M AN US In this paper, we propose a diffusive prey-predator system with mutually interfering predator (Crowley-Martin functional response) and prey reserve. In particular, we develop and analyze both spatially homogeneous model based on ordinary differential equations and reaction-diffusion model. We mainly investigate the global existence and boundedness of positive solution, stability properties of homogeneous steady state, non-existence of non-constant positive steady state, conditions for Turing instability and Hopf bifurcation of the diffusive system analytically. Conventional stability properties of the non-spatial counterpart of the system are also studied. The analysis ensures that the prey reserve leaves stabilizing effect on the stability of temporal system. The prey reserve and diffusive parameters leave parallel impact on the stability of the spatio-temporal system. Furthermore, we illustrate the spatial patterns via numerical simulations, which show that the model dynamics exhibits diffusion controlled pattern formation by different interesting patterns. Introduction CE 1 PT Keywords: Crowley-Martin functional response; Prey reserve; Hopf bifurcation; Global stability; Positive steady states; Turing pattern AC The interactions of individual organisms with each other and the environment are the most important distinctive features in ecological systems. Predation is an important type of interaction that affects the population dynamics of all species. As a result, one of the main objectives of ecologists is to get insight into prey-predator relationship [1, 2, 3, 4]. The rate of prey consumption by an average predator (functional response) is one of the significant characteristics of the prey-predator relationship. Functional response determines stability and bifurcation scenarios with diverse parameter range of ecological systems. The real functional responses not only show the consistency of data fitting with experimental results but also play an important role for accurate predictions of models. Functional responses can primarily be classified in two categories: (i) prey dependent (linear and Holling type functional responses) [5, 6], (ii) predator dependent (ratio dependent, Beddington-DeAngelis, Crowley-Martin and Hassel-Varley) [7]–[13]. The intuitive 1 ACCEPTED MANUSCRIPT = Y CE dY dT X C(1 − m)Y = X R(1 − ) − , K A1 + B1 (1 − m)X + C1 Y + B1 C1 (1 − m)XY PT dX dT ED M AN US CR IP T and experimental observations [14]–[16] indicate that the feeding rate of per unit consumer decreases because of mutual interference among predators. The assumption that predator’s predation decreases due to high predator density (interference among the predator individuals) even when prey density is high (presence of prey handling or searching of prey by predator individual) was modelled mathematically by Crowley and Martin [17] (hereafter Crowley-Martin functional response and the CM model). For the detailed study of Crowley-Martin functional response and its formulation, one can refer [18]–[22] and references cited therein. The experimental study on prey (Paramecium caudatum) and predator (Didinium nasatum) shows that in a completely homogeneous environment, Didinium destructs all prey and it collapses subsequently. Moreover, due to many reasons (such as, environmental pollution, over predation, unregulated harvesting, etc.) several species have became extinct and many others are at the verge of extinction. This also signifies the persistence (the survival of species for a longer time) of the system [23, 24]. Improving conditions or adaptations [25] of natural habitat, imposing restriction on species harvesting, creating natural reserves, establishment of protected areas are the suitable environmental strategies to save the species from getting extinct. More specifically, some fractions of the prey population can partially be protected by using prey refuges (reserves). Environmental heterogeneity also enhances spatial refuges of prey. In a broader sense, any strategy that decreases predation is refuge. Mite predator-prey interactions often exhibit spatial refugia, which afford the prey some degree of protection from predation and reduce the chance of extinction due to predation [26]. The term reserve zone/refuge has received considerable attention in the dynamics of prey-predator model ([27]–[33]). In particular, spatial prey refuge is one of the most relevant behavioural traits that affects the dynamics of the prey-predator systems [34]. Motivated by the above discussion, in the proposed system, we introduce both predator dependent functional response (Crowley-Martin) as well as prey cover (reserve). We investigate the following prey-predator model: −D+ (1) EC(1 − m)X , A1 + B1 (1 − m)X + C1 Y + B1 C1 (1 − m)XY AC where X(t) and Y (t) represent the prey and predator populations at time t, respectively, with X(0) > 0 and Y (0) > 0. Here, it is assumed that in the absence of predator population, prey population follows logistic growth rate while predators are subjected to natural mortality. The system (1) is incorporated by a refuge which protects mX of prey, where m ∈ [0, 1). Thus only (1 − m)X of the prey is available to predator. The parameters C, E, B1 , C1 are also termed as saturating Crowley-Martin type response function parameters. The parameters E and C stand for the effect for capture rate and conversion factor denoting the newly born predators for each captured prey and effect of handling time for predators, respectively. The parameter C1 measures the magnitude of interference among predators. The biological meaning of all other parameters can be found in [35]. All the parameters, R, K, C, D, E, A1 , B1 , C1 in the system (1) assume only positive values and will be considered as constants throughout our discussion. After non-dimensionalization of the system (1) with the following set of variables and parameters: 2 ACCEPTED MANUSCRIPT C B1 K C1 D CEK B1 C1 K X , y = Y, a = ,b = ,c = ,d = ,e = and f = , K RA1 A1 A1 R RA1 A1 the system (1) takes the following form: dx a(1 − m)y = x (1 − x) − = xf1 (x, y), dt 1 + b(1 − m)x + cy + f (1 − m)xy (2) dy e(1 − m)x = y −d+ = yf2 (x, y). dt 1 + b(1 − m)x + cy + f (1 − m)xy CR IP T t = RT, x = PT ED M AN US Living status of population in their natural environment are very rarely spatially homogeneous [36]. Often this heterogeneity takes the form of a patchiness or other type of structural complexity when areas with high population density alternate with areas where the given species is virtually absent depending upon a particular environmental condition [37]. For instance, an insect population living in corn fields in considerable densities may be almost absent in the space between the fields simply because of unfavourable conditions there. In spatially inhomogeneous case, the existence of a non-constant time-independent positive solution (stationary pattern) ensures the dynamical richness of the system. Moreover, local interactions in space can give rise to large scale spatio-temporal patterns (e.g., spiral waves, spatio-temporal chaos (turbulence), stationary (Turing-type) patterns and transitions between these modes). These phenomena are found in many organizational levels of biotic systems, and their occurrence and properties are largely independent of the precise interaction structure. Such type of spatio-temporal patterns have significant consequences for fundamental biological processes. Pattern formation [38]–[43] study in reaction-diffusion systems [44, 45] has been an informative and active research area since the pioneering works of Alan Turing [46] in 1952. Turing patterns of various diffusive prey-predator models have been studied extensively from several researchers [47, 48, 49]. For detailed study of pattern of diffusive prey-predator systems, one can refer [46, 50] and references cited therein. Taking care of random movement of the individuals in two dimensional habitat, we incorporate diffusion terms into the growth equations for prey and predator populations given by (2), to obtain the following reaction diffusion equations: AC CE a(1 − m)y ∂u = u (1 − u) − + d1 ∇2 u, ∂t 1 + b(1 − m)u + cv + f (1 − m)uv (3) ∂v e(1 − m)u = v −d+ + d2 ∇2 v, ∂t 1 + b(1 − m)u + cv + f (1 − m)uv for (X, t) ∈ Ω × R, where Ω ⊂ R2 , and subjected to initial conditions u(X, 0) = u0 (X) ≥ 0, v(X, 0) = v0 (X) ≥ 0, X ∈ Ω (4) and no-flux boundary condition ∂v ∂u = = 0, X ∈ ∂Ω, t > 0. ∂n ∂n (5) Here, we should note that the boundary of the region Ω is in C ∞ i.e., ∂Ω ∈ C ∞ The boundary is smooth enough and there is no break or sharp edge in the boundary. In Eq. 3 ACCEPTED MANUSCRIPT CR IP T (3), the positive parameters d1 and d2 are the diffusion coefficients of prey and predator, ∂2 ∂2 respectively, ∇2 = + is the usual two dimensional Laplacian operator in variable ∂x2 ∂y 2 X ≡ (x, y) ∈ Ω ⊂ R2 and Ω is a square bounded (simply connected) domain with smooth ∂ boundary ∂Ω. n is the outward unit normal vector of the boundary ∂Ω while stands ∂n for outward drawn normal derivative on the boundary ∂Ω. Thus the prey and predator density at location ‘X’ and time ‘t’ is denoted by u(X, t) and v(X, t) respectively. The zero flux boundary conditions (5) assure that there is no fluxes of populations through the boundary, i.e., no external input is imposed from outside. In the present paper, we will be concerned with the following two important questions: i How the prey reserve parameter (m) and diffusion coefficients (d1 , d2 ) affect the dynamical behaviour (co-existence, local and global stability and pattern structure) of the systems (2) and 3). AN US ii To explore the similarities/differences between sufficient parametric conditions for co-existence of steady states, local and global stability of interior equilibrium of temporal system and its spatio-temporal counterpart. 2.1 Analysis of the system (2) Positivity and boundedness of the solutions PT 2 ED M The rest of the paper is structured as follows. Permanence, persistence, local and global stability of the equilibria and related dynamical behaviors of the system (2) are considered in the next section. The dynamics of the associated spatio-temporal model (3) is discussed in the Section 3. More precisely, in Section 3, we study the existence and uniqueness of the global positive solution and nonexistence of nonconstant positive steady states and Turing pattern formation in the system (3). In Section 4, we perform some carefully designed numerical simulations to substantiate our analytical findings. Paper ends with a section on discussion and the ecological significance of our analytical results. Some future directions are also mentioned. CE n Let R+ denotes the set of all non-negative real numbers and R+ = {x ∈ Rn : x = (x1 , ..., xn ) where xi ∈ R+ , ∀ i = 1, 2, ...n}. Here we discuss the invariance of the boundary 2 and the interior of R+ along with conservation of overall energy or biomass flow. AC Lemma 2.1. The positive quadrant (Int(R2+ )) is invariant for the system (2). Proof. The proof is straightforward and interested readers can refer to [51]. Lemma 2.2. All solutions of the system (2) that start in R2+ are confined to the region 1 D = {(x, y) ∈ R2+ : x(t) ≤ 1, 0 ≤ η(t) ≤ 1 + } as t → ∞ for all positive initial value 4d a 2 (x(0), y(0)) ∈ R+ where η(t) = x + y(t). e Remark 2.3. The proof of this theorem follows similarly as in [52]. Geometrically, the Lemma 2.2 assures that all the trajectories of the system (2) initiating from any point in R2+ ultimately lie in the fixed bounded region defined by D. Hence, the flow/dynamical system associated with system (2) is dissipative. 4 ACCEPTED MANUSCRIPT 2.2 Permanence Theorem 2.4. Define M1 := 1, M2 := (e−bd)(1−m), m1 := Any positive solution {x(t), y(t)} of system (2)satisfies m1 ≤ lim inf t→∞ x(t) ≤ lim supt→∞ x(t) ≤ M1 , m2 ≤ lim inf t→∞ y(t) ≤ lim supt→∞ y(t) ≤ M2 , o n d c . whenever m < min 1 − , 1 − a (e − bd)m1 c−a(1−m) , m2 c := m1 (e−bd)(1−m)−d . (c+f m1 )d Remark 2.5. For the detailed proof of the Theorem 2.4, the reader is referred to [53]. CR IP T Using the definitions of lim sup and lim inf, one can claim that there exist four positive constants ti (i = 1, 2, 3, 4) such that m1 ≤ x(t), ∀t ≥ t1 , x(t) ≤ M1 , ∀t ≥ t2 , m2 ≤ y(t), ∀t ≥ t3 , y(t) ≤ M2 , ∀t ≥ t4 . If we define m∗ = min(m1 , m2 ), M ∗ = max(M1 , M2 ) and T = max(ti , 1 ≤ i ≤ 4, i ∈ N), then we have the following permanence result for (2). AN US Theorem 2.6. Suppose that the time independent coefficients a, b, c, d, e, f are positive and satisfy n o c d m < min 1 − , 1 − , (6) a (e − bd)m1 then the model system (2) is permanent. Theorem 2.7. If e < bd, then limt→∞ y(t) = 0. 2.3 ED M Remark 2.8. From the predator equation in (2), the Theorem 2.7 follows immediately. One can easily observe that when f = m = 0 (absence of prey reserve and prey handling at high prey density), the above discussion also remains valid. In this case, the system (2) reduces to the system with Beddington-DeAngelis functional response by Cantrell and Cosner [54] and Takeuchi et al. [55] with density dependence. The theorem 2.4 and 2.6 reduces to the corresponding results in [55]. Existence and stability around biomass equilibria PT The co-existence conditions of all trivial and non-trivial equilibria of the system (2) and their dynamical behaviors are presented in the following table. Feasibility condition (i) E 0 (0, 0) always (ii) E 1 (1, 0) always AC CE Equilibrium and coordinate (iii) E ∗ (x∗ , y ∗ ) Jacobian matrix and eigenvalues 1 0 0 −d λ1 = 1, λ2 = −d a(1−m) −1 − 1+b(1−m) 0 a(1 − m) > (c + f (1 − m)x∗ )(1 − x∗ ) Stability status saddle ! e(1−m) 1+b(1−m) e(1−m) −d + 1+b(1−m) −d + λ1 = −1, λ2 = j11 j12 , j21 j22 λ1 + λ2 = j11 + j22 λ1 λ2 = j11 j22 − j12 j21 stable if m < 1 − d e−bd stable if j11 < 0 Here, x∗ is given by the root of the following cubic equation in z αz 3 + βz 2 + γz + δ = 0, (7) 2 2 where, α = ef (1 − m) , β = e(1 − m) c − f (1 − m) + df (1 − m), γ = a(1 − m) (e − bd) − c d + e(1 − m) , δ = −ad(1 − m). 5 ACCEPTED MANUSCRIPT CR IP T The constant term (δ) of the equation (7) is negative and the leading coefficient (α) is always positive. Thus, the Descartes rule of sign ensures that the above cubic equation (7) possesses at least one positive root, say x∗+ . Furthermore, for this value of x∗+ , the (1 + b(1 − m)x∗+ )(1 − x∗+ ) ∗ ∗ = . will be given by y+ corresponding value of y+ a(1 − m) − (c + f (1 − m)x∗+ )(1 − x∗+ ) a(1 − m)x∗ y ∗ (b(1 − m) + f (1 − m)y ∗ ) , Moreover, j11 = x∗ f1 x = −x∗ + (1 + b(1 − m)x∗ + cy ∗ + f (1 − m)x∗ y ∗ )2 −a(1 − m)x∗ (1 + b(1 − m)x∗ ) j12 = x∗ f1 y = < 0, (1 + b(1 − m)x∗ + cy ∗ + f (1 − m)x∗ y ∗ )2 e(1 − m)y ∗ (1 + cy ∗ ) j21 = y ∗ f2 x = > 0, (1 + b(1 − m)x∗ + cy ∗ + f (1 − m)x∗ y ∗ )2 e(1 − m)x∗ y ∗ (c + f (1 − m)x∗ ) j22 = y ∗ f2 y = − < 0. (1 + b(1 − m)x∗ + cy ∗ + f (1 − m)x∗ y ∗ )2 The following theorem provides alternative conditions for local asymptotic stability of E ∗. AN US Theorem 2.9. If permanence condition (Theorem 2.6) and m < 1 + hold, then the positive equilibrium E ∗ is locally asymptotically stable. f − a + x∗ (bc − f ) ab(c − x∗ ) Global stability: We now discuss the global stability of interior equilibrium solution E ∗ = (x∗ , y ∗ ) by constructing a suitable Liapunov functional. M Theorem 2.10. Let ξ = (1 + b(1 − m)x∗ + cy ∗ + f (1 − m)x∗ y ∗ )(1 + b(1 − m)x + cy + f (1 − m)xy). If permanence conditions given in Theorems (2.4 and 2.6) and m > ay ∗ 1− hold then E ∗ = (x∗ , y ∗ ) is globally asymptotically stable. (1 − x∗ )(b + aM2 ) ED Proof. The system (2) can be re-written in the following form: PT dx a(1 − m)y ∗ = x − (x − x∗ ) + dt 1 + b(1 − m)x∗ + cy ∗ + f (1 − m)x∗ y ∗ a(1 − m)y , − 1 + b(1 − m)x + cy + f (1 − m)xy (8) dy e(1 − m)x e(1 − m)x = y − + . dt 1 + b(1 − m)y ∗ + cy ∗ + f (1 − m)x∗ y ∗ 1 + b(1 − m)y + cy + f (1 − m)xy CE ∗ AC Now, define V (x, y) : R2+ → R, such that. V (x, y) = V1 (x) + V2 (y), (9) where, V1 (x) = x − x∗ − x∗ ln(x/x∗ ) , V2 (y) = η y − y ∗ − y ∗ ln(y/y ∗ ) . The constant η is a positive constant and is defined in the discussion below. The time derivatives of V1 6 ACCEPTED MANUSCRIPT and V2 along the solution of (8) are CR IP T dV1 dx = (x − x∗ ) = (x − x∗ ) − (x − x∗ ) + dt xdt a(1 − m)y ∗ a(1 − m)y , − 1 + b(1 − m)x∗ + cy ∗ + f (1 − m)x∗ y ∗ 1 + b(1 − m)x + cy + f (1 − m)xy dV2 dy = η(y − y ∗ ) = η(y − y ∗ ) − dt ydt e(1 − m)x∗ e(1 − m)x . + 1 + b(1 − m)x∗ + cy ∗ + f (1 − m)x∗ y ∗ 1 + b(1 − m)x + cy + f (1 − m)xy From the time derivative of (9) and a straightforward calculation, we obtain dV dt a(1 − m)2 y ∗ (b + ay) ηe(1 − m)x∗ (c + f (1 − m)x) (x − x∗ )2 − (y − y ∗ )2 ξ(x, y) ξ(x, y) ηe(1 − m)f (1 + cy ∗ ) − a(1 − m)(1 + a(1 − m)x∗ ) + (x − x∗ )(y − y ∗ ). ξ(x, y) = −(x − x∗ )2 + AN US With the help of Theorems (2.4) and (2.6) and choosing η in such a way that ηe(1 − m)(1 + cy ∗ ) = a(1 − m)(1 + (1 − m)ax∗ ), we obtain h h dV b + aM2 ) i (c + f (1 − m)m1 ) i < − 1 − a(1 − m)2 (x − x∗ )2 − ηe(1 − m)x∗ (y − y ∗ )2 .(10) dt ξ(x, y) ξ(x, y) Using the value of ξ(x, y), coefficient of (x − x∗ )2 takes the following form: M a(1 − m)2 (b + aM2 ) (1 + b(1 − m)x∗ + cy ∗ + f (1 − m)x∗ y ∗ )(1 + b(1 − m)x + cy + f (1 − m)xy) (1 − m)(b + aM2 )(1 − x∗ ) , ≤ −1 + (a1 + c1 L) ay ∗ ED − 1+ ay ∗ . Thus we conclude that if all the conditions (1 − x∗ )(b + aM2 ) dV of Theorem (2.10) are satisfied, then < 0 along all the trajectories in R2+ except E ∗ . dt Therefore, E ∗ is globally asymptotically stable. Spatiotemporal model and dynamics AC 3 CE PT which is negative if m > 1− In this section, to capture the dynamics and inhomogeneous distribution of predator and prey, we discuss the diffusive model (3). 3.1 Global existence and boundedness of a positive solution of (3) Definition 3.1. A vector of the form α = (α1 , α2 , . . . , αn ) where each component αi is a nonnegative integer, is called a multilinear of order |α| = α1 + α2 + . . . + αn . For any ∂ |α| α multiindex α, one can define, D U (x) = . ∂xα1 1 , ∂xα2 2 , . . . ∂xαnn 7 ACCEPTED MANUSCRIPT Definition 3.2. Sobolev Space: Let 1 ≤ p < ∞. The Sobolev space W k,p (U ) consists of locally summable functions u : U → R such that for each multiindex α with |α| ≤ k, Dα u exists in weak sense and belongs to Lp (U ). Here if p = 2 use notation H k (U ) = W k,2 (U ) (k = 0, 1, 2...). The system of coupled partial differential equations (3) can be written as ∂w = G(w) + M ∆w, ∂t (11) CR IP T where w = (u, v)T and M = diag{d1 , d2 }, G(w) = (G1 (u, v), G2 (u, v)), G1 (u, v) = u (1 − a(1−m)y e(1−m)u u) − 1+b(1−m)u+cv+f (1−m)uv , G2 (u, v) = v − d + 1+b(1−m)u+cv+f (1−m)uv . Moreover, let H 2 (Ω) denotes Sobolev space of all square integrable functions such that all the derivatives up to second order are also in L2 (Ω). Here the initial conditions u(X, 0) = u0 (X), and v(X, 0) = v0 (X), belongs to Sobolev space H 2 (Ω) ∪ 0. We define Aw = −M ∆w for w ∈ D(A), where A(w) = {w : w ∈ H 2 (Ω)×H 2 (Ω)}. Hence using this notation, the system (11) can be rewritten in the following abstract differential equation: AN US ∂w + Aw = G(w). (12) ∂t One can easily establish the existence of a local solution for (12) and the non-negativity of components u and v of the solution of initial-boundary value problem (3). The proof of uniform boundedness of the solution of (3) suffices the global existence of the positive solution. For this, we recall the following lemma: M Lemma 3.3. Let w(X, t) satisfies the following equation PT ED ∂w − ∆w = w(1 − u) X ∈ Ω, ∂t ∂w = 0 X ∈ ∂Ω, t > 0, ∂n u(X, 0) = u0 (X) > 0, X ∈ Ω. Then lim u(X, t) = 1 for each X ∈ Ω. t→∞ AC CE Theorem 3.4. (Weak Persistence) Let Ω̄ = Ω ∪ ∂Ω. The non-negative solution (u, v) satisfy lim sup max u(., t) ≤ 1, lim sup max v(., t) ≤ t→∞ u∈Ω̄ t→∞ u∈Ω̄ (1 − m)(e − bd) . d (13) Proof. Using positivity of u and v, we obtain ∂u = d1 ∆u + u(1 − u). ∂t Applying Lemma 3.3, ∃ a T > 0 and an arbitrary small > 0 such that u(X, t) < 1 + in Ω × [T, ∞). Similarly, from second equation of (3), after using the bound of u, one can also estimate the bound for v lim sup max v(., t) ≤ t→∞ v∈Ω̄ e(1 − m)(1 + ) − bd(1 − m) . cd 8 ACCEPTED MANUSCRIPT Moreover, the arbitrariness of gives that lim sup max v(., t) ≤ v∈Ω̄ t→∞ (1 − m)(e − bd) . cd The above discussion can also be summarized in the following theorem: Theorem 3.5. The system (3) has a unique, nonnegative and bounded solution w = (u, v)T such that CR IP T w ∈ C([0, ∞), X) ∩ C 1 ([0, ∞), X) ∩ C([0, ∞), D(A)). Theorem 3.6. (Strong Persistence) The system (3) is strongly persistent i.e. lim inf max u(., t) > 0 and lim inf max v(., t) > 0 if t→∞ t→∞ u∈Ω̄ v∈Ω̄ (1 − m)(e − bd)(c − a(1 − m)) > d, o n (1−m)(e−bd)(c−a(1−m))−d . where 0 = min c−a(1−m) , c c+f (c−a(1−m))d AN US c > a(1 − m), Proof. Proof follows similarly as in Theorem 3.4. 3.2 M Remark 3.7. Theorems 3.4 and 3.6 also ensure that conditions and bounds for the persistence of the temporal system (system (1) and Theorem 2.6) are same as results obtained for the system with diffusion (system (3) Theorem 3.4). Homogeneous steady state and its stability analysis ED One can easily establish the existence and uniqueness of a positive homogeneous steady state U ∗ = (u∗ , v ∗ ) of the system (3). To investigate the local stability of U ∗ and the effect of diffusion on the system (3), we linearize the system (3) around U ∗ . Substituting T PT u(X, t), v(X, t) = U ∗ + Υ(X, t) = U ∗ + T U (X, t), V (X, t) , (14) AC CE in the system (3) and picking up all the terms which are linear in U and V, the system (3) takes the following from: ∂Υ = D∆Υ + J(U ∗ )Υ, ∂t where D = diag(d1 , d2 ). More precisely, ∂U = j11 U + j12 V + d1 ∆U, ∂t ∂V = j21 U + j22 V + d1 ∆V, ∂t (15) where jij for i, j = 1, 2 are same as obtained for the model system without diffusion (refer 2 the section 2.3). Let λi ∞ i=1 denotes the set eigenvalues of −∆ on Ω ∈ R with zero flux boundary condition such that 0 < λ1 < λ2 < · · · . Further, let χ(λi ) denote the eigenspace corresponding to λi , and ψij , for j = 1, 2, 3, · · · , dimχ(λi ), be an orthonormal basis of 9 ACCEPTED MANUSCRIPT ∂v ∂u χ(λi ). Let Hij = {cψij ; c ∈ R2 } and Y = {(u, v) ∈ C 1 (Ω̄) × C 1 (Ω̄) : ∂n = ∂n = 0}. Ldimχ(λi ) L∞ Consider the following direct sum Xi = i=1 Hij and Y = i=1 Xi . Then the space Xi is invariant under the operator P = D∆ + J(U ∗ ) for each i ≥ 1. Moreover, λ is an eigenvalue of P if and only if it is also an eigenvalue of −d1 λi + j11 j12 ∗ ∗ Ai = J(λi , U ) = −λi D + J(U ) = . j21 −d2 λi + j22 The eigenvalues of Ai are determined by where (16) CR IP T η 2 − Bη + C = 0, B = −λi (d1 + d2 ) + j11 + j22 = −λi (d1 + d2 ) + trJ(U ∗ ), C = d1 d1 λ2i − (d1 j22 + d2 j11 )λi + detJ(U ∗ ). (17) (18) AN US It may be noted here that J(U ∗ ) is same as JE ∗ (refer the Jacobian matrix associated with temporal system). From equations (16)-(18) and using Routh-Hurwitz criteria [56], one can summarize the above discussion in the following theorem: Theorem 3.8. (i) The positive equilibrium U ∗ of the system (3) is locally asymptotically if B > 0 and C > 0. M (ii) If j11 < 0 (refer the subsection 2.3), then the positive equilibrium is locally asymptotically stable in the presence as well as in the absence of diffusion. ED Remark 3.9. The sufficient condition (Section 2.3) for the asymptotic stability of positive equilibrium E ∗ of the temporal system (2) also implies the asymptotic stability of homogeneous steady state solution U ∗ of the diffusive system (3). In the next theorem, we discuss the global stability of positive constant steady-sate U of (3). The existence of globally stable constant steady-state signifies two important points: (i) Predator and prey species can maintain steady state for all future time under certain parametric restriction. (ii) It also ensures the non-existence of spatial patterns. PT ∗ AC CE Theorem 3.10. The positive homogeneous steady state U ∗ of the system (3) is globally asymptotically stable if the corresponding equilibrium solution of (2) is globally asymptotically stable. Proof. To prove the Theorem 3.10, we consider the following Lyapunov functional Z V (t) = W1 (u(X, t), v(X, t))dX, where Ω u W1 (u, v) = u − u − u log ∗ + K u ∗ ∗ Z v v∗ η − v∗ dη. η Here K is a positive constant and is defined in discussion below. Differentiating W1 with ∗ v−v ∗ 1 1 = u−u , ∂W = ve(1−m) . The time derivative of V (t) respect to u and v, we obtain ∂W ∂u u ∂v 10 ACCEPTED MANUSCRIPT CR IP T along the trajectory of the system (3) gives Z ∂W1 ∂u ∂W1 ∂v dV = + dX dt ∂u ∂t ∂v ∂t Ω Z i u − u∗ h a(1 − m)uv = d1 ∆u + u(1 − u) − u 1 + b(1 − m)u + cv + f (1 − m)uv Ω i ∗ h v−v e(1 − m)uv + d2 ∆v − dv + dX e(1 − m)v 1 + b(1 − m)u + cv + f (1 − m)uv Z v − v∗ u − u∗ d1 ∆u + d2 ∆v = u e(1 − m)v Ω Z h i e(1 − m)v ∗ (u − u ) (1 − u) − + 1 + b(1 − m)u + cv + f (1 − m)uv Ω i ∗ h v−v e(1 − m)u + −d+ dX (19) e(1 − m) 1 + b(1 − m)u + cv + f (1 − m)uv = I1 (t) + I2 (t), AN US where I1 (t) and I2 (t) denote first and second integral in (19), respectively. Using homogeneous boundary conditions, we obtain Z u − u∗ v − v∗ I1 (t) = d1 ∆u + d2 ∆v dX u e(1 − m)v Ω Z d2 v ∗ d1 u∗ 2 2 |∇u| + |∇v| dX ≤ 0. =− (20) u2 cdv 2 Ω ED M Further, after some straightforward calculations and following the steps involved in the proof of the Theorem 2.10, one can also show that I2 (t) < 0. This completes the proof of Theorem 3.10. Nonexistence of nonconstant positive steady state CE 3.3 PT Remark 3.11. The Theorem 3.10 also concludes that if all the conditions of the Theorem 2.10 are well satisfied, then the homogeneous steady state of the associated spatiotemporal model (3) is globally asymptotically stable. AC For the non-existence of non-constant positive equilibrium of the system (3), we establish the non-existence of non-constant positive equilibrium of the following second order elliptic system: −d1 ∆u = u (1 − u) − −d2 ∆v = v − d + a(1 − m)v , 1 + b(1 − m)u + cv + f (1 − m)uv e(1 − m)u 1 + b(1 − m)u + cv + f (1 − m)uv with X ∈ Ω and zero-flux boundary condition ∂u ∂v = = 0, X ∈ ∂Ω. ∂n ∂n 11 , (21) ACCEPTED MANUSCRIPT Before proving main theorem for non-existence of non-constant positive steady state, we introduce two important lemmas and one theorem related to priori lower and upper bounds for the positive solution of (21). Lemma 3.12. (Harnack Inequality) [57] Let m(X) ∈ C(Ω). Further let ϕ(X) ∈ C 2 (Ω) ∩ C 1 (Ω) be a positive solution of ∆ϕ(X) + m(X)ϕ(X) = 0 in Ω subject to ∂ϕ = 0 on ∂Ω. ∂ν Then ∃ a positive constant M∗ = M∗ (km(X)k∞ ) such that max ϕ(X) ≤ M∗ min ϕ(X). Ω Ω (a) If π ∈ C 2 (Ω) ∩ C 1 (Ω) satisfies ∆π(X) + ϕ(X, π(X) ≥ 0 in Ω, and π(X0 ) = max π(X), then ϕ(X0 , π(X0 )) ≥ 0. Ω ∂π ≤ 0 on ∂Ω, ∂ν AN US (b) If π(X) ∈ C 2 (Ω) ∩ C 1 (Ω) satisfies CR IP T Lemma 3.13. [57] Let ϕ ∈ C(Ω × R). ∆π(X) + ϕ(X, π(X) ≤ 0 in Ω, ∂π ≥ 0 on ∂Ω, ∂ν (22) (23) and π(X0 ) = min π(X), then ϕ(X0 , π(X0 )) ≤ 0. Ω M Remark 3.14. The standard regularity theory of elliptic system [61] ensures that the positive solution of (21) is contained in C 2 (Ω) × C 2 (Ω). Hence one can use the above lemmas to find out some proper estimates (lower and upper bounds) for the solution of (21). (i) Any positive solution (u, v) of (21) satisfy ED Theorem 3.15. max u(X) ≤ 1, Ω (1 − m)(e − bd) , cd (24) PT Ω max v(X) ≤ CE (ii) Assume that d∗ be any fixed positive constant such that d1 , d2 ≥ d∗ and e(1 − m) > d(1 + b(1 − m)). Then ∃ a positive constant Q = Q(a, b, d, d∗ , e, m) such that min u(X) ≥ Q, Ω min v(X) ≥ Q. (25) Ω AC Proof. Both the inequalities of (24) are followed immediately by using Lemma 3.13 and equation (21). Furthermore, we denote u(X1 ) = min u(X), u(X2 ) = max u(X), (26) v(Y1 ) = min v(Y ), v(Y2 ) = max v(Y ). (27) Ω Ω Ω Ω 12 ACCEPTED MANUSCRIPT Using the Lemma 3.13, equations (26) and (27) in (21), we obtain the following inequalities: a(1 − m)v(X1 ) , 1 + b(1 − m)u(X1 ) + cv(X1 ) + f (1 − m)u(X1 )v(X1 ) e(1 − m) − d(1 + b(1 − m)) u(Y1 ) v(Y1 ) ≥ , d(1 + f (1 − m)u(Y1 )) e(1 − m) − d(1 + b(1 − m)) u(Y2 ) , v(Y2 ) ≤ d(1 + f (1 − m)u(Y2 )) 1 − u(X1 ) ≤ (28) (29) (30) AN US CR IP T where e(1 − m) > d(1 + b(1 − m)). Equations (29) and (30) together with (26) and (27) imply that e(1 − m) − d(1 + b(1 − m)) u(Y1 ) v(Y1 ) ≥ d(1 + f (1 − m)u(Y1 )) e(1 − m) − d(1 + b(1 − m)) u(X1 ) ≥ , (31) d(1 + f (1 − m)u(X1 )) e(1 − m) − d(1 + b(1 − m)) u(Y2 ) v(Y2 ) ≤ d(1 + f (1 − m)u(Y2 )) e(1 − m) − d(1 + b(1 − m)) u(X2 ) ≤ . (32) d(1 + f (1 − m)u(X2 )) Using (27) and (32) in (28), we obtain M 1 ≤ u(X1 ) + M ∗ u(X2 ), (33) e(1 − m) − d(1 + b(1 − m)) . d h i a(1−m)v −1 Now we define m(X) = d1 1 − u − 1+a(1−m)u+cv+f (1−m)uv . Then applying the Harnack inequality (Lemma 3.12) on ∆u(X) + m(X)u(X) = 0 (from first equation of (21)), we have (34) max u(X) ≤ Q1 min u(X), i.e., u(X2 ) ≤ Q1 u(X1 ), PT ED where M ∗ = a(1 − m) X∈Ω X∈Ω AC CE where Q1 is a positive constant which depends on km(X)k∞ . Combining (34) with (33), we have min u(X) = u(X1 ) ≥ (1 + M ∗ Q1 )−1 := Q2 . (35) X∈Ω From (31), we have min v(y) = v(y1 ) ≥ y∈Ω M ∗ (1 + M ∗ Q)−1 := Q3 . a(1 − m)(1 + f (1 − m)(1 + M ∗ Q)−1 ) (36) Defining Q = min{Q2 , Q3 }, we have min u(X) = u(X1 ) ≥ Q, min v(Y ) = v(Y1 ) ≥ Q. Theorem 3.16. Suppose that d2 ∗ > e(1−m)(1+b(1−m))−d be a fixed positive constant. Then λ1 one can obtain a positive constant d1 ∗ such that the model system (21) has no non-constant positive solution provided d1 > d1 ∗ and d2 > d2 ∗ . 13 ACCEPTED MANUSCRIPT Proof. Let us denote F1 (u, v) = u(1 − u) − a(1−m)uv , 1+b(1−m)u+cv+f (1−m)uv e(1−m)uv . 1+b(1−m)u+cv+f (1−m)uv R 1 ũ = Ω Ω u(X)dX, ṽ F2 (u, v) = −dv + Let (u, v) be a positive solution of (21). Moreover, we define R = Ω1 Ω v(X)dX, for any u, v ∈ L1 (Ω) and κ(u, v) = (1 + b(1 − m)u + cv + f (1 − m)uv)(1 + b(1 − m)ũ + cṽ + f (1 − m)ũṽ). Integrating the first equation of (21) after multiplying by (u − ũ), we obtain Z 2 F1 (u, v)(u − ũ)dX = {F1 (u, v) − F1 (ũ, ṽ)}dX Ω ZΩ h uv = (u − ũ) − (u2 − ũ2 ) − a(1 − m) 1 + b(1 − m)u + cv + f (1 − m)uv Ω i ũṽ dX − 1 + b(1 − m)ũ + cṽ + f (1 − m)ũṽ Z h (v + cvṽ)(u − ũ)2 ≤ (u − ũ)2 − a(1 − m) κ(u, v) Ω (ũ + b(1 − m)uũ)(u − ũ)(v − ṽ) i + dX. (37) κ(u, v) d1 |∇(u − ũ)| dX = CR IP T Ω Z Z AN US Using the bounds of u and v from (25) and P1 = ab in (37), we have Z Z h i 2 d1 |∇(u − ũ)| dX ≤ (u − ũ)2 + 2P1 |u − ũ||v − ṽ| dX. Ω (38) Ω Similarly integrating the second equation of (21) after multiplying by (v − ṽ), we obtain Z Z F2 (u, v)(u − ũ)dX = {F2 (u, v) − F2 (ũ, ṽ)}dX Ω Z h uv = − d(v − ṽ) + e(1 − m) 1 + b(1 − m)u + cv + f (1 − m)uv Ω i ũṽ − dX. (39) 1 + b(1 − m)ũ + cṽ + f (1 − m)ũṽ d2 |∇(v − ṽ)| dX = Ω PT ED Ω 2 M Z 2 AC CE Now by defining P2 = 1 + c(1 − m)2 (e−bd) and using the bounds of u and v from (24) in c2 d 2 (39), we obtain Z Z 2 d2 |∇(v − ṽ)| dX ≤ (e(1 − m)(1 + b(1 − m)) − d)(v − ṽ)2 Ω Ω +e(1 − m)P2 (u − ũ)(v − ṽ) dX. (40) Let us define R = P2 + 2P1 . Now from inequalities (38) and (40), we have Z Z Z 2 2 {d1 |∇(u − ũ)| + d2 |∇(v − ṽ)| }dX ≤ (u − ũ)2 + (e(1 − m)(1 + b(1 − m)) Ω Ω 2 −d)(v − ṽ) + 2R|u − ũ||v − ṽ| dX. (41) 14 ACCEPTED MANUSCRIPT Now, using the Poincare inequality, the above inequality (41) yields Z Z Z R 2 2 d1 λ1 |∇(u − ũ)| + d2 λ1 |∇(v − ṽ)| dX ≤ 1+ (u − ũ)2 + e(1 − m)(1 δ Ω Ω +b(1 − m)) − d + δR (v − ṽ)2 dX, (42) CR IP T for an arbitrary positive constant δ. As d2 λ1 > e(1 − m)(1 + b(1 − m)) − d, one can find a sufficiently small δ ∗ such that d2 λ1 ≥ e(1 − m)(1 + b(1 − m)) − d + δ ∗ R. Moreover, if we ∗ define d1 = λ1 1 + δR∗ , then d1 λ1 − 1 + δR∗ > 0 (refer the condition of Theorem 3.16). Hence one can conclude that u = ũ and v = ṽ. 3.4 AN US Remark 3.17. The Theorem 3.16 determines threshold values of diffusivity constants for the non-existence of a non-constant positive solution of the system (21). The condition d2 ∗ > e(1−m)(1+b(1−m))−d also ensures that the prey reserve parameter (m) affects the nonλ1 existence of a non-constant positive solution which in turn implies that the prey refuge may leave some effect for pattern formation. Turing space and Hopf instabilities As we know that any solution of the linear system (15) can be expanded into a Fourier series. Hence γ1 mn (t) sin k(km , kn )X, m,n=0 ∞ X m,n=0 ∞ X ∞ X γ2 mn (t) sin k(km , kn )X, ∞ X unm (X, t) = ED M U (X, t) = V (X, t) = m,n=0 vnm (X, t) = (43) m,n=0 AC CE PT where km and kn are corresponding wavenumbers. Now, replacing U and V from unm vnm into (15), wee have dγ1 nm = (j11 − d1 k 2 )γ1 nm + j12 γ2 nm , dt (44) dγ2 nm = j21 γ1 nm + (j22 − d2 k 2 )γ2 nm , dt 2 where k 2 = kn2 + km . Equation (44) is a system of two linear differential equations and a general solution of (44) has the form A exp(λ1 t) + B exp(λ2 t). Here the constants A and B are determined with the help of initial conditions (4) while λ1 and λ2 are the eigenvalues of the associated coefficient matrix J= j11 − d1 k 2 j12 j21 j22 − d2 k 2 15 . ACCEPTED MANUSCRIPT Here λ1 and λ2 are the solutions of the characteristic equation associated with the matrix J λ2 − tr(J) + det(J), (45) where tr(J) = tr(JE ∗ )−(d1 +d2 )k 2 and det(J) = det(D)k 4 +(−d2 j11 −d1 j22 )k 2 +det(JE ∗ ). It is obvious that the diffusive instability holds, if tr(J) < 0 and det(J) < 0 and j11 < 0 are satisfied. Moreover, if j11 < 0 holds, then tr(JE ∗ ) < 0 and hence tr(J) < 0. The condition det(J) < 0 is sufficient for the diffusive instability and which is equivalent to the inequality CR IP T S(k 2 ) = det(D)k 4 + (−d2 j11 − d1 j22 )k 2 + det(JE ∗ ) < 0. (46) Here S is a quadratic function with respect to k 2 . Moreover the minimum value of S(k 2 ) is obtained at 1 2 d1 j22 + d2 j11 . (47) k 2 = kmin = 2 det(D) AN US 2 Hence more stronger condition for diffusive instability becomes S(kmin ) < 0 i.e., 1 d2 j11 + d1 j22 > 2(det(D)det(JE ∗ ) 2 . (48) As a consequence, we have the following theorem. Theorem 3.18. The system (3) undergoes diffusion-induced instability if the condition of Theorem 2.9, relations (46) and (48) hold, i.e., 1 (49) Numerical simulations ED 4 M d2 j11 + d1 j22 > 2(det(D)det(JE ∗ ) 2 . PT In order to show the analytical findings and stability results obtained in the previous sections graphically, we numerically simulate the solutions of both the model systems (diffusive and non-diffusive) with different sets of parameter values. AC CE Case 1: The Non-diffusive System: Here we execute two interesting aspects, namely coexistence and dynamics of the system (2) with respect to two different significant parameters, viz., prey refuge (m) and maximum predation rate (a). The parameter values for the system (2) are listed in the Table 1. First fixing the parameter a = 20 and other parameters from Table 1, when we vary prey refuge parameter m throughout its whole range [0, 1) (see bifurcation diagram with respect to m: Fig.1(i)), then we observe that at the lower range of m coexistence equilibrium point (E ∗ ) of the system (2) is unstable focus and eventually system converges to nearest stable limit cycle around E ∗ (Fig.1(iα): phase-plane diagram with m = 0.1). At the intermediate range of m, system (2) converges to its coexistence equilibrium point (E ∗ ) (Fig.1(iβ): phase-plane diagram with m = 0.6). Here, the system (2) changes its unstable focus to stable focus status via Hopf-bifurcation [58, 59, 60] with respect to parameter m. At the higher range of m, due to the lack of food (prey), predator goes to extinction while prey goes to its highest density (carrying capacity) by gradually decreasing foraging efficiency of predator (Fig.1(iγ: phase-plane diagram with m = 0.97). 16 ACCEPTED MANUSCRIPT CR IP T Now coexistence equilibrium point (E ∗ ) becomes infeasible (ecological point of view) and predator-free equilibrium point (E 1 ) becomes stable node. In this case, a transcritical bifurcation occurs between E ∗ and E 1 [62]. On the other hand, keeping fixed the parameter m = 0.01 and other parameters from Table 1, if we vary maximum predation rate parameter a throughout a large range [2, 20] (refer the bifurcation diagram with respect to a: Fig.1(ii)), then we notice that at the lower range of a, system (2) converges to its coexistence equilibrium point (E ∗ ) (Fig.1(iiα): phase-plane diagram with a = 4). Moreover, at higher range of a, coexistence equilibrium point (E ∗ ) of the system (2) becomes an unstable focus and eventually system converges to nearest stable limit cycle around E ∗ (Fig.1(iiβ): phase-plane diagram with a = 15). Here, system (2) changes its stable focus to unstable focus status via Hopf-bifurcation with respect to parameter a. Case 2: The Diffusive System: M AN US For the diffusive system, we will open up discussion from three different perspectives: (1) suppose that there is no prey refuge (i.e., m = 0) and observe the effect of prey diffusion on the pattern formation; (2) assume that there is some prey refuge (i.e., m 6= 0) and examine the effect of prey diffusion on the pattern formation; (3) observation on the system with different intensities of prey refuge along with other system parameters remain fixed. In every situation, the relation between S(k 2 ) and k is shown and the corresponding pattern structure are also demonstrated. In the process of numerical simulation, we select the Neumann boundary condition with a size of 200 × 200. Generally speaking, Neumann boundary is also called as zeroflux boundary, which implies that the domain is isolated or insulated from the external environment. We assume that the space step is dx = 1 and the time step is dt = 0.005. The initial distributions of both prey and predator are taken as follows: v(x, y, 0) = v ∗ + ξ2 , ED u(x, y, 0) = u∗ + ξ1 , (50) CE PT where ξi (i = 1, 2) is a tiny random perturbation term. We run the simulations until pattern reaches a stationary state or does not seem to change its main characteristics any longer. Simultaneously, we fix d2 = 50 and other parameter values are given in Table 1 except the parameters a, m, d1 . Here, we just show the pattern formation of prey population (u) in the xy-plane since predator population (v) exhibits extremely similar distributions. AC Case 2.1: Without refuge and the effect of diffusion on pattern structure In this subsection, we suppose that the prey have no any refuges (m = 0) and observe the effect of prey diffusion on the pattern formation. From the equation (46), it is obvious that the condition for the diffusion-induced instability is equivalent to S(k 2 ) < 0. According to the Fig. 2(i), it is clear that the interval of negativity of the function S(k 2 ) becomes more smaller as d1 increases i.e., the corresponding diffusive system will be unstable only when the prey diffusion coefficient is smaller and Turing pattern emerges under only this circumstance. The pattern evolution with different prey diffusion coefficients are shown in Figure 3. When the value of d1 is smaller, the center of region is occupied by shortstripe pattern and surrounded by hot-spot pattern with high density (Fig. 3(i)). If d1 takes a middle value, the density of center is uniform and the uniform density region is 17 ACCEPTED MANUSCRIPT surrounded by fewer short-stripe pattern and the outermost part of the area is hot-spot pattern with intermediate density (Fig. 3(ii)). When the value of d1 is high, the pattern structure is same as the above situation, but the uniform density central region is more larger and the hot-spot pattern becomes more sparse (Fig.3(iii)). It is obvious that the diffusion coefficient of the prey reduces the range of pattern region, that is to say, the diffusion coefficient have a stabilizing effect on the spatiotemporal system (3). Case 2.2: With refuge and the effect of diffusion on pattern structure PT ED M AN US CR IP T In this case, we consider a constant prey refuge (m = 0.1) and examine the effect of prey diffusion on the pattern structure. The plot of S(k 2 ) as a function is shown in Fig.2(ii), which is similar to Fig.2(i). The only difference is that both the interval of negativity of the function S(k 2 ) and the corresponding functional value becomes more smaller. Namely, when the prey refuge is taken into account, prey diffusion can still give rise to Turing instability, however the prey refuge seemingly reduces the instability probability. The associated pattern formation is shown in Fig.4. When d1 is smaller, one can observe that the density of central zone is uniform and surrounded by hot-spot pattern (Fig.4(i)). If we increase d1 slightly, analogous-spiral pattern emerges and it is encompassed by cold-spot pattern (Fig.4(ii)). When d1 is higher, the cold-spot pattern fades away and the region is occupied by analogous-spiral pattern (Fig.4(iii)). Moreover, the central region appears the most safest region for the prey and the farther from the center, the lower the security seems. From the perspective of the preys and predators, in order to avoid being attacked by predators, the preys migrate towards the center. Moreover, when the preys are protected at a certain percentage m, (1 − m) of preys is only available to predator. The protective effect of population aggregation will lose a certain superiority and the preys have to escape at a higher speed to avoid being eaten. This ensures that the diffusion coefficient of preys will be higher because Turing pattern is induced by diffusion, or more strictly speaking, it is caused by the value of dd12 . When the refuge is considered, the value of dd21 becomes smaller, to some extent, which reduces the possibility of Turing pattern. It is in complete agreement with the result of the Fig.2(ii). Case 2.3: The effect of the refuge on pattern structure AC CE In this situation, we fix the diffusion coefficient of prey at d1 = 0.01 and vary m from 0.01 to 0.1 and then we examine the effect of prey refuge on the pattern structure. The plot of S(k 2 ) with different values of refuge m is shown in Fig. 2(iii). The Fig. 2(iii) demonstrates that the interval of negativity of the function S(k 2 ) becomes more smaller as m increases. Furthermore, if the value of m is greater than 0.14, the functional value of S(k 2 ) will be positive and the Turing instability does not occur. Comparing Fig. 2(iii) with Fig. 2(i), we can conclude that the prey refuge and diffusion have the parallel impact on the diffusion-induced instability. The stationary-state patterns for different values of m are shown in Fig. 5 and the pattern structure is almost consistent to Fig.3. This perfectly illustrates that prey refuge and diffusion coefficient have the parallel impact on pattern formation and structure. Thus one can summarize the above discussion by saying that prey reserve and diffusive coefficients both leave an stabilizing effect on the spatiotemporal system (3). 18 ACCEPTED MANUSCRIPT 5 Discussion AC CE PT ED M AN US CR IP T Many efforts have been put to manage sub-population of the vulnerable cheetas (Acinonyx jubatus) in small reserves. Like vulnerable cheetas there are many other species those are either driven to extinction or at the verge of extinction, due to overexploitation, over predation, indiscriminate harvesting, environmental pollution etc. There are several biological examples where individual prey can decrease their predation risk by using a refuge, which may be habitat that is spatially separate, a phase of prey life history that is temporarily separate, a collective behaviour of prey such as aggregation into herds. This motivates us for the creation of biospheres reserve zones/refuges and this may decrease the interaction among species. While creating protected areas for a species several factors should be taken into consideration like it’s capacity to house number of individuals of the species to be protected, dynamics of the ecosystems supporting these species, economics associated with the maintenance of the protected area, bio-economics of the ecosystem and several others. The establishment of protected areas (reserve zones) is the potential solution for optimal harvesting problems. In this paper, we have investigated the dynamical behaviors of a prey-predator system (both for diffusive and non-diffusive) with Crowley-Martin functional response and prey reserve. The prey refuge is also a significant factor as far as the size of habitat and co-existence of species are concerned. The temporal system is thoroughly examined with different types of bifurcations. Parametric conditions of co-existence and local asymptotic stability of each biomass equilibrium are presented systematically. We have derived threshold values of prey reserve parameter (m) for the permanence and global stability of the temporal system (2) (refer the Theorems 2.6 and 2.10). The parameter ranges of existence of multiple bifurcations are identified. More precisely, bifurcation diagram with respect to prey reserve (m) and maximum predation rate (a) ensures the occurrence of transcritical bifurcation and Hopf-bifurcation (refer the Figure 1). A rigorous investigation of the global dynamics and bifurcations of patterned solutions of the spatio-temporal counterpart of the prey-predator system under zero-flux boundary conditions is also presented. We have pointed out the similarity/difference between the sufficient parametric conditions obtained for persistence, local and global stability of both the temporal and its diffusive counterpart. The theoretical analysis reveals the fact that global stability of the interior equilibrium solution (E ∗ ) of non-diffusive system (2) is sufficient for the global asymptotic stability of positive homogeneous state (U ∗ ) of diffusive system (3). The role of prey reserve in the non-existence of non-constant solution has been investigated which suggests that prey reserve must leave some effect on pattern formation (see the Theorem 3.16). We have carried out extensive numerical simulations to explore the effect of prey reserve (m) and diffusion coefficients on pattern structure. The obtained results ensure that the prey reserve and diffusion coefficient leave parallel impact on pattern structure. In particular, the presence of prey reserve reduces the possibility of Turing pattern. In brief, the value of the present study can be summarized in three folds. First: it analyzes the effects of prey refuge and maximum predation rate parameters on the stability of the non-diffusive system. System shows diverse dynamical phenomena with respect to these two significant system parameters and system switches its stability via Hopfbifurcation and changes stability with other feasible system equilibrium via transcritical bifurcation. Second: the spatial system classifies diffusion-driven instability conditions and gives corresponding theorems. Third: it illustrates spatial patterns close to the 19 ACCEPTED MANUSCRIPT Acknowledgments AC 6 CE PT ED M AN US CR IP T onset of instability bifurcation via numerical simulations, indicating that the dynamics of the diffusive system exhibits a prey refuge and diffusion controlled pattern formation. Moreover, as our study shows that the effect of prey refuge for pattern formation in two dimensional spatial domain has profound effect on predator-prey interaction. On one side, this may enrich the spatial dynamics of the effect of prey refuge on the predatorprey model while on the other side it may provide important information to practical ecologists about the impact of prey refuge on pattern formation. The stabilizing role of prey reserve in the dynamics and stability of prey-predator system have been discussed by few researchers [52, 53] but not much with Crowley-Martin functional response [54]. As the present study deals with combined effects of diffusion and prey refuge, it also generalizes the existing works in [52, 53, 54]. The results obtained for global stability and permanence in [52, 53], can easily be obtained from our study just by ignoring mutual interference among predator species. Future directions: It should be noted that pattern selection is an important method to study reaction-diffusion equation, which has played a great role in chemistry, spatial epidemiology, and ecosystems [63, 64, 65]. Predator-prey relationship is very crucial in population and abundant pattern structures can easily be found in predator-prey systems, such as spot, stripe, mixed, labyrinth, etc [65, 66]. In the biological systems, the pattern formation has a very significant effect on the stability of community and ecosystem. For instance, the spot pattern indicates that the prey (or predator) exhibits patchy distribution and population aggregation (the preys avoid predation and the predators add the chance of catching preys). In this sense, for the pattern formation in the predator-prey system, the two significant issues for further investigations are: (1) previous studies in predator-prey system have focused on the local role of population, however, the nonlocal effects of population have a great impact on themselves for the population movement. As a result, the population model with nonlocal effect should be more reasonable [67]. The travelling wave solutions and entire solutions of the nonlocal models have received some attention, but no one has ever done the research of pattern dynamics, which will be a new perspective for population dynamics. (2) In order to maintain the ecosystem stability, ecologists suggest several suitable measures, such as provision of prey refuge, supply of additional food and so on. We wonder that what kind of control measures are most reasonable or profitable. Maybe, we can combine the pattern formation with control theory of the reaction-diffusion equation to find the optimal control strategy and the corresponding optimal pattern structures. The project is funded by the National Natural Science Foundation of China under Grant 11671241, 131 Talents of Shanxi University, Program for the Outstanding Innovative Teams (OIT) of Higher Learning Institutions of Shanxi, and Natural Science Foundation of Shanxi Province Grant No. 201601D021002. References [1] A.A. Berryman, The origin and evolution of predator-prey theory, Ecology 73(5) (1992) 1530-1535. 20 ACCEPTED MANUSCRIPT [2] C. Xu, Y. Shao, Bifurcations in a predator-prey model with discrete and distributed time delay, Nonlinear Dynamics 67(3) (2012) 2207-2223. [3] G.-Q. Sun, C.-H. Wang, L.-L. Chang, Y.-P. Wu, L. Li, Z. Jin, Effects of feedback regulation on vegetation patterns in semi-arid environments, Appl. Math. Model. 61 (2018) 200-215. [4] C. Xu, X. Tang, M. Liao, M., Stability and bifurcation analysis of a delayed predatorprey model of prey dispersal in two-patch environments, Appl. Math. Comput 216(10) (2010) 2920-2936. CR IP T [5] C.S. Holling, Some characteristics of simple types of predation and parasitism, Canada Entomology 91 (1959) 385-398. [6] T.K. Kar, H. Matsuda, Global dynamics and controllability of a harvested preypredator system with Holling type III functional response, Nonlinear Analysis: Hybrid Systems 1 (2007) 59-67. AN US [7] S. Liu, E. Beretta, A stage-structured predator-prey model of Beddington-DeAngelis type , SIAM J. Appl. Math. 66(4) (2006) 1101-1129. [8] J.R. Beddington, Mutual interference between parasites or predators and it’s effect on searching efficiency, J. Anim. Ecol. 44(1) (1975) 331-340. [9] S.B. Hsu, T.W. Hwang, Y. Kuang, Global analysis of the Michaelis-Menten-type ratio-dependent predator-prey system, J. Math. Biol. 42 (2001) 489-506. M [10] C. Xu, P. Li, Oscillations for a delayed predator-prey model with Hassell-Varley-type functional response, Comptes Rendus Biologies 338(4) (2015) 227-240. ED [11] J.P. Tripathi, S.S. Meghwani, M. Thakur, S. Abbas, A modified Leslie-Gower predator-prey interaction model and parameter identifiability, Commun. Nonlinear Sci. Numer. Simulat. 54 (2018) 331-346. PT [12] S. Ruan, D. Xiao, Global analysis in a predator-prey system with non-monotonic functional response, SIAM J. Appl. Math. 61(4) (2001) 1445-1472. AC CE [13] J.P. Tripathi, S. Tyagi, S. Abbas, Global analysis of a delayed density dependent predator-prey model with Crowley-Martin functional response, Commun. Nonlin. Sci. Numer. Simulat. 30 (2016) 45-69. [14] M.P. Hassell, Mutual interference between searching insect parasites, J. Anim. Ecol. 40 (1971) 473-486. [15] D. Jana, J.P. Tripathi, Impact of generalist type sexually reproductive top predator interference on the dynamics of a food chain model, In. J. Dynam. Control (2016) DOI: 10.1007/s40437-016-0255-9. [16] B. Zimmermann, H.K. Sand, P. Wabakken, O. Liberg, H.P. Andreassen, Predatordependent functional response in wolves: from food limitation to surplus killing, J. Anim. Ecol. 84 (2015) 102-112. 21 ACCEPTED MANUSCRIPT [17] P.H. Crowley, E.K. Martin, Functional responses and interference within and between year classes of a dragonfly population, J. N. Amer. Benth. Soc., 8 (1989) 211-221. [18] G.T. Sklaski, J.F. Gilliam, Functional responses with predator interference: viable alternative to Holling type II model, Ecology 82 (11) (2001) 3083-3092. [19] R.D. Parshad, A. Basheer, D. Jana, J.P. Tripathi, Do prey handling predators really matter: Subtle effects of a Crowley-Martin functional response, Chaos, Solitons and Fractals 103 (2017) 410-421. CR IP T [20] R.K. Upadhyay, S.N. Raw, V. Rai., Dynamic complexities in a tri-trophic food chain model with Holling type II and Crowley-Martin functional response, Nonlin. Anal.: Model. Control 15(3) (2010) 361-375. [21] J.P. Tripathi, Almost periodic solution and global attractivity for a density dependent predator-prey system with mutual interference and Crowley-Martin response function, Differ. Equ. Dyn. Syst. (2016) DOI: 10.1007/s12591-016-0298-6. AN US [22] H. Yin, X. Xiao, X. Wen, K. Liu, Pattern analysis of a modified Leslie-Gower predator-prey model with Crowley-Martin functional response and diffusion, Comput. Math. Appl. 67 (2014) 1607-1621. [23] V. Hutson, K. Schmitt, Permanence and the dynamics of biological systems, Math. Biosci. 111 (1992) 1-71. M [24] S. Abbas, M. Banerjee, N. Hungerbuhler, Existence, uniqueness and stability analysis of allelopathic stimulatory phytoplankton model, J. Math. Anal. Appl. 367 (2010) 249-259. ED [25] V. Krivan, D. Jana, Effects of animal dispersal on harvesting with protected areas, J. Theor. Biol. 364 (2015) 131-138. PT [26] X. Guan, W. Wang, Y. Cai, Spatiotemporal dynamics of a Leslie-Gower predatorprey model incorporating a prey refuge, Nonlinear Anal. RWA 12 (2011) 2385-2395. [27] J. Smith, Models in ecology, Cambridge University Press, Cambridge, 1974. AC CE [28] D. Jana, R. Agrawal, R.K. Upadhyay, Dynamics of generalist predator in a stochastic environment: effect of delayed growth and prey refuge, Appl. Math. Comput. 268 (2015) 1072-1094. [29] D. Jana, S. Ray, Impact of physical and behavioral prey refuge on the stability and bifurcation of Gause type Filippov prey-predator system. Model. Earth Syst. Environ. 2(1) (2016) 2-24. [30] B. Dubey, A prey-predator model with a reserved area, Nonlinear Analysis: Modelling and Control 12(4) (2007) 479-494. [31] B. Dubey, P. Chandra, P. Sinha, A model for fishery resource with reserve area, Nonlinear Analysis. RWA. 4 (2003) 625-637. 22 ACCEPTED MANUSCRIPT [32] Y. Huang, F. Chen, L. Zhong, Stability analysis of a predator-prey model with Holling type III response function incorporating a prey refuge, Appl. Math. Comput. 182 (2006) 672-683. [33] R. Yang, J. Wei, Stability and bifurcation analysis of a diffusive predator-prey system in Holling type III with a prey refuge, Nonlinear Dynamics 79 (2015) 631-646. [34] D. Jana, Chaotic dynamics of a discrete predator-prey system with prey refuge, Appl. Math. Comput. 24 (2013) 848-865. CR IP T [35] D.L. DeAngelis, R.A. Goldstein, R.V. O’Neill, A model for tropic interaction, Ecological Society of America 56(4) (1975) 881-892. [36] L.A.D. Rodrigues, D.C. Mistro, S. Petrovskii, Pattern formation, long-term transients, and the Turing-Hopf bifurcation in a space-and time-discrete predator-prey system, Bull. Math. Biol. 73(8) (2011) 1812-1840. AN US [37] R.S. Cantrell, C. Cosner, Spatial ecology via reaction-diffusion equations, Wiley, West Sussex, 2003. [38] M. Banarjee, S. Petrovskii, Self-organized spatial pattern and chaos in a ratio dependent predator-prey system, Theoretical Ecology 4 (2011) 37-53. [39] G.-Q. Sun, A. Chakraborty, Q.X. Liu, Z. Jin, K.E. Anderson, B.L. Li, Influence of time delay and nonlinear diffusion on herbivore outbreak, Commun. Nonlinear Sci. Numer. Simulat. 19(5) (2014) 1507-1518. M [40] R.K. Upadhyay, V. Volpert, N.K. Thakur, Propagation of Turing patterns in plankton model, J. Biol. Dynam. 6 (2012) 524-538. ED [41] Y. Wang, J. Cao, G.-Q. Sun, J. Li, Effect of time delay on pattern dynamics in a spatial epidemic model, Physica A 412 (2014) 137-148. PT [42] X-P. Yan, C-H Zhang, Stability and turing instability in a diffusive predator-prey system with Beddington-DeAngelis functional response, Nonlinear Analysis: RWA. 20 (2014) 1-13. CE [43] X.C. Zhang, G.-Q., Sun, Z. Jin, Spatial dynamics in a predator-prey model with Beddington-DeAngelis functional response, Physical Review E 85(2) (2012) 021924. AC [44] H. Merdan, S. Kayan, Hopf bifurcations in Lengyel–Epstein reaction–diffusion model with discrete time delay, Nonlinear Dynam. 79(3) (2015) 1757-1770. [45] W. Zuo, Y. Song, Stability and bifurcation analysis of a reaction–diffusion equation with distributed delay, Nonlinear Dynamics 79(1) (2015) 437-454. [46] A.M. Turing, The chemical basis of mokphogenesis, Philos. Trans. R. Soc. Lond. 237 (1952) 37-72. [47] M. Banerjee, S. Abbas, Existence and Non-existence of spatial patterns in a ratiodependent predator-prey model, Ecol. Complex. 21 (2015) 199-214. 23 ACCEPTED MANUSCRIPT [48] M. Banerjee, S. Banerjee, Turing instabilities and spatio-temporal chaos in ratiodependent Holling-Tanner model, Math. Biosci. 236 (2012) 64-76. [49] G. Zhang, W. Wang, X. Wang, Coexistence states for a diffusive one-prey and twopredators with B-D functional response, Comput. Math. Appl. 387 (2012) 931-948. [50] L. Xue, Pattern formation in a predator-prey model with spatial effect, Physica A 391 (2012) 5987-5996. CR IP T [51] J.P. Tripathi, S. Abbas, M. Thakur, Local and global stability analysis of two prey one predator model with help, Commun. Nonlin. Sci. Numer. Simulat. 19 (2014) 3284-3297. [52] L. Li, C.-H. Wang, S.-F. Wang, M.-T. Li, L. Yakob, B. Cazelles, Z. Jin, W.-Y. Zhang, Hemorrhagic fever with renal syndrome in China: Mechanisms on two distinct annual peaks and control measures, Internat. J. Biomath. 11 (2018) 1850030. AN US [53] J.P. Tripathi, S. Abbas, M. Thakur, Dynamical analysis of a prey–predator model with Beddington-DeAngelis type function response incorporating a prey refuge, Nonlinear Dynam. 80 (2015) 177-196. [54] R.S. Cantrell, C. Cosner, On the dynamics of predator-prey models with the Beddington-DeAngelis functional response, J. Math. Anal. Appl. 257 (2001) 206-222. [55] H. Li, Y. Takeuchi, Dynamics of the density dependent predator-prey system with Beddinton-DeAngelis functional response, J. Math. Anal. Appl. 374 (2011) 644-654. M [56] C. Xu, Y. Wu, Bifurcation and control of chaos in a chemical system, Appl. Math. Model. 39(8) (2015) 2295-2310 ED [57] Y. Lou, W-M. Ni, Diffusion self-diffusion and cross diffusion, J. Differ Equations 131 (1996) 79-131. PT [58] C. Xu, Stability and bifurcation analysis in a viral model with delay, Math. Meth. Appl. Sci. 36(10) (2013) 1310-1320. CE [59] R. Curtu, Singular Hopf bifurcations and mixed-mode oscillations in a two-cell inhibitory neural network, Physica D 239(9) (2010) 504-514. AC [60] C. Xu, X. Tang, M. Liao, Frequency domain analysis for bifurcation in a simplified tri-neuron BAM network model with two delays, Neural Networks 23(7) (2010) 872880. [61] D. Gilbarg, N.S. Trudinger, Elliptic partial differential equations of second order, Springer-Verlag, New York, 1996. [62] D. Jana, N. Bairagi, Habitat complexity, dispersal and metapopulations: Macroscopic study of a predator–prey system, Ecol. Complex. 17 (2014) 131-139. [63] Q. Ouyang, Patterns formation in reaction diffusion systems. Shanghai Sci-Tech. Education Publishing House, Shanghai, 2000. 24 ACCEPTED MANUSCRIPT [64] G.-Q. Sun, C.-H. Wang, Z.-Y. Wu, Pattern dynamics of a Gierer–Meinhardt model with spatial effects, Nonlinear Dynamics 88 (2017) 1385-1396. [65] G.-Q. Sun, Z.-Y. Wu, Z. Wang, Z. Jin, Influence of isolation degree of spatial patterns on persistence of populations, Nonlinear Dynam. 83 (2016) 811-819. [66] A. Yochelis, Y. Tintut, L.L. Demer, A. Garfinkel, The formation of labyrinths, spots and stripe patterns in a biochemical approach to cardiovascular calcification. New J. Phys. 10 (2008) 055002. CR IP T [67] X.-Q. Zhao, J. Borwein, P. Borwein, Dynamical systems in population biology, Springer, New York, 2003. AN US Parameters Symbols Prey saturation constant b Predation interference constant c Predator mortality rate d Conversion factor e Predator interference (in the presence of prey handling) f Default values 2.10 1.10 0.2 5 0.001 AC CE PT ED M Table 1: System parameters and their default values chosen for simulation 25 ACCEPTED MANUSCRIPT (i) a=20 1 0.9 0.9 (iα) m=0.1 y 0.8 0.8 0 0.7 0 0 0.8 0.7 (iβ) m=0.6 0 0.4 (iiβ) a=15 1 y y 1 0.6 0.6 0.5 0 0 1 CR IP T 0 0.5 (iγ) m=0.97 1 y Population (iiα) a=4 1 y 1 (ii) m=0.01 1 x 0.4 0 x 0 1 0.4 x x 0 1 0.3 0.2 0.2 y AN US 0.3 y 0.1 0 0.1 0 0.2 0.4 0.6 0.8 1 0 2 5 10 15 20 a M m AC CE PT ED Figure 1: Bifurcation diagram of system (2) w.r.t m (Fig.(1)(i)) and three phase-plane diagrams w.r.t lower (Fig.(1)(iα)), intermediate (Fig.(1)(iβ)) and higher (Fig.(1)(iγ)) value of m keeping a = 20. Bifurcation diagram of system (2) w.r.t a Fig.(1)(ii) and three phase-plane diagrams w.r.t lower (Fig.(1)(iiα)) and higher (Fig.(1)(iiβ)) value of a keeping m = 0.01. Other parametric values are given in Table 1. (i) (ii) (iii) Figure 2: The graph of the function S(k 2 ) for parametric values given in Table 1 and for (i) m = 0, a = 5, d2 = 50, (ii) m = 0.1, a = 5, d2 = 50 and (iii) a = 5, d1 = 0.01, d2 = 50. 26 ACCEPTED MANUSCRIPT (i) (ii) (iii) (ii) (i) AN US CR IP T Figure 3: Snapshots of contour pictures of the time evolution of prey (u) in the xy-plane with d2 = 50 at 300000 iterations where (i): d1 = 0.01; (ii): d1 = 0.02; (iii): d1 = 0.03 and other parametric values are given in Table 1. (iii) (i) CE PT ED M Figure 4: Snapshots of contour pictures of the time evolution of prey (u) in the xy-plane with m = 0.1, d2 = 50 at 300000 iterations where (i): d1 = 0.012; (ii): d1 = 0.015; (iii): d1 = 0.018 and other parametric values are given in Table 1. (ii) (iii) AC Figure 5: Snapshots of contour pictures of the time evolution of prey (u) in the xy-plane with d1 = 0.01 and d2 = 50 at 300000 iterations where (i): m = 0.01; (ii): m = 0.05; (iii): m = 0.1 and other parametric values are given in Table 1. 27

1/--страниц