Explosive tritrophic food chain models with interference: A comparative study

Explosive tritrophic food chain models with interference: A comparative study

Journal Pre-proof

Explosive tritrophic food chain models with interference: A comparative study Debaldev Jana, Ranjit Kumar Upadhyay, Rashmi Agrawal, Rana D. Parshad, Aladeen Basheer PII: DOI: Reference:

S0016-0032(19)30853-1 https://doi.org/10.1016/j.jfranklin.2019.11.049 FI 4284

To appear in:

Journal of the Franklin Institute

Received date: Revised date: Accepted date:

28 February 2018 16 September 2019 26 November 2019

Please cite this article as: Debaldev Jana, Ranjit Kumar Upadhyay, Rashmi Agrawal, Rana D. Parshad, Aladeen Basheer, Explosive tritrophic food chain models with interference: A comparative study, Journal of the Franklin Institute (2019), doi: https://doi.org/10.1016/j.jfranklin.2019.11.049

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. 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. © 2019 Published by Elsevier Ltd on behalf of The Franklin Institute.

Highlights • The effect of top predator interference on the dynamics of a tritrophic food chain model is the principal goal to study general framework for three species food chain models in which intermediate and top predators are specialist and generalist types respectively where the top predator grows by sexual reproduction. • This concept is used to design five different models in this paper by five different combinations of functional responses of specialist intermediate predator (food uptake process follows either prey dependent, Holling type III/IV or prey-predator dependent, BD-functional response) and sexually reproductive generalist top predator (food uptake follows strictly prey-predator dependent, BD/CM functional response) respectively assuming the fact that the growth of intermediate predator is mediated by gestation delay. • Stability of each equilibria of the non-delayed counterpart of each model and the existence of Hopf-bifurcation for the coexisting equilibrium point of all delayed systems are established. • The delayed and non-delayed models can blow-up in finite time under sufficient conditions on the initial data. Furthermore, some dynamic behaviors of the systems are clarified by some numerical tests. The numerical simulations show that the gestation delay can act as a damping mechanism and prevent blow-up in certain critical range of gestation period. • Ecological implications of this phenomenon are discussed.

1

Explosive tritrophic food chain models with interference: A comparative study Debaldev Jana1 , Ranjit Kumar Upadhyay2 , Rashmi Agrawal3 , Rana D. Parshad4 and Aladeen Basheer5 1 Department

of Mathematics & SRM Research Institute,

SRM Institute of Science and Technology, Kattankulathur - 603 203, Tamil Nadu, India 2 Department

of Mathematics & Computing,

Indian Institute of Technology (Indian School of Mines), Dhanbad - 826 004, Jharkhand, India 3 Department

of Science & Humanities,

Indian Institute of Information Technology Dharwad, Hubli-580029, Karnataka, India 4 Department

of Mathematics,

lowa State University, Ames, IA 50011, USA 5 Department

of Mathematics,

University of Georgia, Athens, GA 30602, USA Abstract Depending upon the choice of food, availability of resource and growth structure, food uptake process of higher trophic level species are significantly complicated and gives interesting dynamical impacts on community food chain. The effect of top predator interference on the dynamics of a tritrophic food chain model is the principal goal to study general framework for three species food chain models in which intermediate and top predators are specialist and generalist types respectively where the top predator grows by sexual reproduction. This concept is used to design five different models in this paper by five different combinations of functional responses of specialist intermediate predator (food uptake process follows either prey dependent, Holling type III/IV or prey-predator dependent, Beddington-DeAngelis (BD)-functional response) and sexually reproductive generalist top predator (food uptake follows strictly prey-predator dependent, BD/Crowley-Martin (CM) functional response) respectively assuming the fact that the growth of intermediate predator is mediated by gestation delay. We establish the stability of each equilibrium point of the non-delayed counterpart of each model and investigate the existence of Hopf-bifurcation for the coexistence equilibrium point of all delayed systems. We also show that both the delayed and non delayed models can blow-up in finite time under sufficient conditions on the initial data. Furthermore, some dynamic behaviors of the systems are clarified by some numerical tests. The numerical simulations show that the gestation delay can act as a damping mechanism and prevent blow-up in certain critical range of gestation period. Ecological implications of this phenomenon are discussed. Keywords : Holling type III and IV functional responses, Beddington-DeAngelis (BD) functional response, Crowley-Martin (CM) functional response, specialist predator, generalist predator, sexual reproduction, finite time blow-up, damping.

2

1

Introduction

Ecological models are increasingly used to emphasize the linkage and relationship between social and physical environments. They demonstrate our responses to such factors over entire life cycles. They are easy to understand and interpret, they are often based on causality and they can include time delay, age/stage structure, functional response, etc. The effect of these factors has illustrated complicated and rich dynamics such as stability, bifurcation, periodic solutions, persistence, etc [1]. Research results herein are helpful to predict the developing/decaying tendency of specific populations to determine the key factors behind extinction and to seek optimum strategies of preventing and controlling extinction of various species [2]. In population dynamics, a functional response refers to the change in the density of prey consumed per unit time by per predator, as the prey density changes [3]. Based on different backgrounds and experimental data, various forms of functional responses have been developed mimicing various processes of energy transfer in predator-prey systems [4]. Categories which specify the functional response: prey dependent (function of only prey density, viz. Holling type I - IV [5]) and prey-predator dependent (function of both prey-predator density, viz. Beddington-DeAngelis [6], Crowley-Martin [7], Hassel-Varley [8] etc) those are describing the predator’s interference on their hunting procedure. The inherent feature of these functions means that the more prey in the environment the better off the predator, which is true in many predator-prey interactions [9]. Also, non-monotonic response occurs at the microbial level when the nutrient concentration reaches a high level an inhibitory effect on the specific growth rate may occur, Holling type IV functional response are used to model such inhibitory effect. Some researchers suggest that a more suitable general predator-prey model should have a prey-predator dependent functional response [10, 11, 12, 13] and all the prey-predator dependent functional responses can provide better description of predator feeding over a range of predator-prey abundances present [14]. Feeding procedure of higher trophic level species is more complicated than its below trophic level species [15, 16]. In general it is shown in the nature that if lower trophic level species is specialist type, then there is high chance to be generalist type of the higher trophic level species [15, 16]. Sexual reproduction is very common process for the growth of that kind of higher trophic level generalist species [17, 18, 19, 20, 21]. Also gestation period is very important factor in population dynamics and assuming that reproduction of predator after consuming prey is not instantaneous but mediated by some time lag required for gestation [22, 23]. It has been shown that many models in the class of where the top predator is modeled via the modified Leslie-Gower scheme for sexually reproductive generalist predator [24, 22]. We first provide some ecological background of invasive species. We are motivated primarily to control non-native species, which is a cental problem in spatial ecology. Data on species such as the invasive Burmese python (Python bivittatus) in the Florida everglades, show an exponential increase in python population, which have resulted in local prey population reducing severely [25]. An invasive species is formally defined as any species capable of propagating itself in a nonnative environment and thus establish a self-sustained population. In the United States alone damages caused by invasive species to agriculture, forests, fisheries and businesses, have been estimated to be $120 billion a year [26]. Another potential issue 3

is that the environment may turn favorable for a certain species while becoming unfavorable for its competitors or natural enemies. This can result in the favored species to outbreak [27]. For example, in the European Alps certain seasonal environmental conditions enable the population of the larch budmoth to become large enough to defoliate entire forests [28]. then for sufficiently large initial data of top predator density, a potential to blow-up/explode in finite time can occur [24]. Biological control is an adopted strategy to limit harmful populations [29]. The objective of a biological control is to establish a management strategy that best controls and decreases the harmful population to healthy levels as opposed to high and risky levels. Naturally, how does one define high level, and further, how well does the biological control actually work, at various high levels? We have recently started investigating this question via the mathematical property of finite time blow-up [30, 31, 32]. This ideas have also been applied to model invasive populations that seems to be ”exploding”, under a variety of ecological scenarios. In the current manuscript we ask: • What is the foraging pressure of intermediate specialist predator on prey if it is modeled by Lotka-Volterra scheme and generalist sexually reproductive top predator on intermediate predator if the top predator is modeled according to the modified Leslie-Gower scheme? • What are the dynamical outcomes when the predation process of intermediate and top predator are considered by different combinations of prey dependent and preypredator dependent functional responses? With respect to the different combinations of functional responses different models will come for analysis. • We ask what would the effects be, if there is predator interference as well. • What will be periodic effect of gestation delay of intermediate predator on the steady state of the system? • In particular, we ask what would the finite time blow up be, if the top generalist predator grows by sexual reproduction as well. • How would the blow-up/blow up-prevention be effected in the case of various functional responses/feeding rates of the preys. In this study, the effect of top predator interference on the dynamics of a general tritrophic food chain model with five different combinations of functional responses is investigated. The manuscript is organized as follows: Models are formulated in Section 2 and stability properties of these models are studied in Section 3. Also, existence of Hopf-bifurcation for coexisting equilibrium point of all models are investigated in this section. In Section 4, finite time blow up behavior is investigated and numerical simulation has been carried out in Section 5. A comparative study between different models and a brief conclusion are in Section 6 and 7 respectively.

4

2

Model formulation

First, we represent prey, intermediate specialist predator and sexually reproductive generalist top predator as X, Y and Z respectively in a general model, the system is:   dX X = a1 X 1 − − F1t .Y, dt K dY (1) = −a2 Y + θF1t−τ .Y − F2 (Y, Z).Z, dt dZ w3 Z 2 = cZ 2 − , dt Y +D where all parameters are positive constants. a1 is the intrinsic growth rate and K is the environmental carrying capacity of X. a2 and θ are the death rate and food conversion efficiency of Y . c denotes the growth rate of Z. D represents the residual loss in Z population due to severe scarcity of its favorite food Y . F1t represents functional response of Y upon X at time t, it may be prey dependent: F1t (X(t)) or prey-predator dependent: F1t (X(t), Y (t)) (when we consider mutual interference of Y ). When we consider the gestation delay (τ ) of Y , then the numerical response becomes θF1t−τ [either F1t−τ (X(t−τ )) or F1t−τ (X(t−τ ), Y (t−τ )) and it is a function at time scale (t − τ )] where θ is the food conversion efficiency of Y . F2 (Y, Z) represents functional response of Z population upon Y and mutual interference of Z is considered here. The ecological basis in deriving the above model is as follows: • Prey (X) grows logistically. • The intermediate specialist predator (Y ) and sexually reproductive generalist top predator (Z) causes a negative feedback on prey (X) and intermediate predator (Y ) [it is clear from predation terms in prey (X) and intermediate predator (Y ) equations −F1t and −F2 (Y, Z)]. • Due to the delayed numerical response, there is a positive feedback term (+θF1t−τ ) in the intermediate predator (Y ) equation at time t, where τ is the gestation time and θ is the food conversion efficiency of Y . • Top predator (Z) is generalist type, so the top predator (Z) is modeled according to the modified LeslieGower scheme. When there is no intermediate predator i.e., Y = 0, then Z has a positive growth rate (c − wD3 ) and can survive in the environment. • In the formulation the functional response of the top predator is not similar to that of its prey (Y ). The generalist top predator grows quadratically as Z 2 , signifying population growth is directly proportional to the product of males and females (Z × Z = Z 2 ) 3 [17, 18, 19, 20]. The population decays due to intraspecies competition as −( Yw+D )Z 2 . Thus if the prey (here Y ) population is large, Z has enough prey, and so the competition coefficient is small. On the other hand if the prey population is small, Z has lack of sufficient prey and so the competition coefficient is large, inducing greater competition amongst the predators Z. The D shows that Z is a true generalist predator and can switch to alternate prey, in case its favorite prey Y goes extinct [21]. 5

Now we consider the feeding abilities of intermediate predator (Y ) upon its prey X and top predator (Z) upon its prey Y that are measured by different functional responses in different five combinations and finally we get five different models in the following forms: (i) Model 1: Combination of Holling III (1st predator) and BD (2nd predator) functional responses Here, we consider the predation of intermediate (Y ) and top predator (Z) follow Holling III and BD functional responses respectively. Therefore, the model system becomes   dX X wX 2 Y , = a1 X 1 − − 2 dt K X +b dY w1 (X(t − τ ))2 Y (t − τ ) w2 Y Z = −a2 Y + − , 2 dt (X(t − τ )) + b Y + βZ + γ w3 Z 2 dZ = cZ 2 − . dt Y +D

(2)

(ii) Model 2: Combination of Holling III (1st predator) and CM (2nd predator) functional responses Here, we consider the predation of intermediate (Y ) and top predator (Z) follow Holling III and CM functional responses respectively. Therefore, the model system becomes   X wX 2 Y dX = a1 X 1 − − 2 , dt K X +b dY w1 (X(t − τ ))2 Y (t − τ ) w2 Y Z = −a2 Y + − , 2 dt (X(t − τ )) + b ρY Z + δY + β1 Z + γ1 dZ w3 Z 2 = cZ 2 − . dt Y +D

(3)

(iii) Model 3: Combination of Holling IV (1st predator) and BD (2nd predator) functional responses Here, we consider the predation of intermediate (Y ) and top predator (Z) follow Holling IV and BD functional responses respectively. Therefore, the model system becomes   dX X = a1 X 1 − − dt K

X2 j

wXY , +X +a

dY w1 X(t − τ )Y (t − τ ) w2 Y Z = −a2 Y + X 2 (t−τ ) − , dt + X(t − τ ) + a Y + βZ + γ j 2

dZ w3 Z = cZ 2 − . dt Y +D

6

(4)

(iv) Model 4: Combination of Holling IV (1st predator) and CM (2nd predator) functional responses Here, we consider the predation of intermediate (Y ) and top predator (Z) follow Holling IV and CM functional responses respectively. Therefore, the model system becomes   dX X = a1 X 1 − − dt K

X2 j

wXY , +X +a

dY w1 X(t − τ )Y (t − τ ) w2 Y Z = −a2 Y + X 2 (t−τ ) − , dt + X(t − τ ) + a ρY Z + δY + β1 Z + γ1

(5)

j 2

dZ w3 Z = cZ 2 − . dt Y +D

(v) Model 5: Combination of BD (1st predator) and CM (2nd predator) functional responses Here, we consider the predation of intermediate (Y ) and top predator (Z) follow BD and CM functional responses respectively. Therefore, the model system becomes   X wXY dX = a1 X 1 − − , dt K X + βY + γ w1 X(t − τ )Y (t − τ ) w2 Y Z dY = −a2 Y + − , dt X(t − τ ) + βY (t − τ ) + γ ρY Z + δY + β1 Z + γ1 dZ w3 Z 2 = cZ 2 − . dt Y +D

(6)

The initial conditions of all the above systems (2)-(6) are given as X(θ) = φ1 (θ) ≥ 0, Y (θ) = φ2 (θ) ≥ 0, Z(θ) = φ3 (θ) ≥ 0, θ ∈ [−τ, 0], φi (0) > 0 (i = 1, 2, 3), where φ : [−τ, 0] → R3 with norm ||φ|| = sup {|φ1 (θ)|, |φ2 (θ)|, |φ3 (θ)|}, such that φ = (φ1 , φ2 , φ3 ). −τ ≤θ≤0

By the fundamental theory of functional differential equations [33], we know that there is a unique solution (X(t), Y (t), Z(t)) to systems (2)-(6) with above initial conditions. In these models, a prey population X serves as the only food for the specialist predator population of Y . This population, in turn, serves as a favourite food for the generalist predator population Z. In this model, K, a, a1 , a2 , j, b, c, d, D, β, γ, ρ, δ, w, w1 , w2 , w3 are positive constants. The interaction between specialist predator Y and the generalist predator Z is modeled by the modified Leslie-Gower scheme where the loss in predator population is proportional to the reciprocal of per capita availability of its most favourite food. We assume that in the absence of the predator the prey population density grows according to a logistic law with carrying capacity K(K > 0) and with an intrinsic growth rate constant a1 (a1 > 0). a2 is the intrinsic death rate of the predator Y in the absence of the only food X. w, w2 , w3 are the maximum values which per capita reduction rate of X, Y, Z can attain respectively, where as w2 is the maximum value which per capita growth rate of Y can attain. a and b are considered as half saturation constant for Holling type III and IV functional response respectively 7

in the absence of any inhibitory effect. The parameter j is the measure of the predator’s immunity from or tolerance of the prey. γ and γ1 represent the protection provided to the prey by its environment for BD and CM functional responses respectively. β and β1 are the intensity of interference between individuals of predator respectively for BD and CM functional responses. δ measures the magnitude of interference among the prey individuals for CM functional response. ρ is considered as the inter species interference between prey and predator for CM functional response. D represent the residual loss in Z population due to severe scarcity of its favourite food Y . c is the growth rate of the generalist predator Z due to sexual reproduction and Z 2 signifies the fact that mating frequency is directly proportional to the number of male as well as female individuals. The mechanism of sexual reproduction of top predator is given in [17, 18, 19, 20]. Now if we rewrite the systems (2)-(6) in general form, then we get: dX = f1 (X, Y, 0), dt dY = f2 (0, Y, Z) + f¯2 (X(t − τ ), Y (t − τ ), 0), dt dZ = f3 (0, Y, Z). dt

(7)

Where all the functions f1 , f2 , f¯2 , f3 are generalized forms of the proposed five models given in Table 1. Model (7) is written in this way only for the analytic analysis purpose which is followed by forthcoming sections.

3

Stability properties of the systems

Feasibility or biologically positivity studies aim to objectively and rationally uncover the strengths and weaknesses of an existing proposed model in a given environment. Ecological stability can refer to types of stability in a continuum ranging from regeneration via resilience (returning quickly to a previous state), to constancy to persistence. The precise definition depends on the ecosystem in question, the variable or variables of interest and the overall context. In the context of conservation ecology, stable population are often defined as ones that do not go extinct. Researchers applying mathematical models from system dynamics usually use Lyapunov stability. Local stability indicates that a system is stable over small short-lived disturbances. Linearize the system (7) at any arbitrary equilibrium point E(x, y, z), we get V˙ (t) = AV (t) + BV (t − τ ) (8)       X(t) f1X f1Y 0 0 0 0 ¯ ¯      Y (t) 0 f2Y f2Z f2X f2Y 0  , where where V (t) = , A = and B = Z(t) 0 f3Y f3Z 0 0 0 all the entries of the Jacobian matrices A and B are given in Table 2. The characteristic equation of (8) is given by λ3 + P λ2 + Qλ + R + e−λτ (Sλ2 + T λ + U ) = 0, 8

(9)

where

3.1

P = −(f1X + f2Y + f3Z ), Q = f2Y f3Z − f2Z f3Y + f1X f2Y + f1X f3Z , R = f1X f2Z f3Y − f1X f2Y f3Z , S = −f¯2Y , T = f3Z f¯2Y − f1Y f¯2X + f1X f¯2X , U = f1Y f3Z f¯2X − f1X f3Z f¯2Y .

Stability properties of the non-delayed counterpart of systems (1-5)

Model systems (2)-(6) have six possible non-negative equilibria, namely E0 (0, 0, 0), EX (K, 0, 0), ˜ EXY (X, ˆ Yˆ , 0), EXZ (K, 0, Z) ¯ and E∗ (X ∗ , Y ∗ , Z ∗ ). The characteristic equation EZ (0, 0, Z), (9) of the system (7) with τ = 0 (non-delayed counterpart of the systems (2)-(6)) at any arbitrary equilibrium point E(x, y, z) is given by: λ3 + (P + S)λ2 + (Q + T )λ + (R + U ) = 0.

(10)

Now we show the feasibility and stability properties of the first five equilibria for the nondelayed (i.e., τ = 0) counterpart of systems (2)-(6) in Table 3-7 along with the corresponding feasibility of coexistence equilibrium point E∗ . Stability condition of the coexistence equilibrium point (E∗ (X ∗ , Y ∗ , Z ∗ )) of the non-delayed counterpart (τ = 0) of the systems (2)-(6) will describe after the ”Stability properties of Model 5”. (i) Stability properties of Model 1 ˜ EXY (X, ˆ Yˆ , 0) Feasibility and stability conditions of E0 (0, 0, 0), EX (K, 0, 0), EZ (0, 0, Z), ¯ of system (2) are given in Table 3. The coexistence equilibrium point of and EXZ (K, 0, Z) Model (2) is E∗ (X ∗ , Y ∗ , Z ∗ ), where Y∗ =

(Y ∗ + γ)[ωKa2 Y ∗ − ω1 a1 X ∗ (K − X ∗ )] ω3 − cD , Z∗ = c ω1 βa1 X ∗ (K − X ∗ ) − ωKY ∗ (a2 β + ω2 )

such that K > X ∗ and ω3 > cD and X ∗ is the positive root of the cubic equation   ωK(ω3 − cD) ∗ ∗3 ∗2 + b X ∗ − Kb = 0. f (X ) = X − KX + ca1 2

3 −cD) > 0 as ω3 > cD. We note that 0 < X ∗ < K. Also, we have f (0) < 0 and f (K) = ωK (ω ca1 Since f (0)f (K) < 0, there is a positive root of this equation lies in (0, K). Therefore, the coexistence equilibrium point E∗ (X ∗ , Y ∗ , Z ∗ ) exist under the following conditions:

X ∗ < K, ω3 > cD.

9

(ii) Stability properties of Model 2 ˜ EXY (X, ˆ Yˆ , 0) Feasibility and stability conditions of E0 (0, 0, 0), EX (K, 0, 0), EZ (0, 0, Z), ¯ and EXZ (K, 0, Z) of system (3) are given in Table 4. The coexistence equilibrium of Model (3) is E∗ (X ∗ , Y ∗ , Z ∗ ), where Y∗ =

(δY ∗ + γ)(X ∗2 (a2 − ω1 ) + a2 b) ω3 − cD , Z∗ = c (ρY ∗ + β)(X ∗2 (ω1 − a2 ) + a2 b) − ω2 (X ∗2 + b)

such that K > X ∗ and ω3 > cD and X ∗ is the positive root of the cubic equation   ωK(ω3 − cD) ∗ ∗3 ∗2 f (X ) = X − KX + + b X ∗ − Kb = 0. ca1 2

3 −cD) We note that 0 < X ∗ < K. Also, we have f (0) < 0 and f (K) = ωK (ω > 0 as ω3 > cD. ca1 Since f (0)f (K) < 0, there is a positive root of this equation lies in (0, K). Therefore, the non-zero equilibrium point E∗ (X ∗ , Y ∗ , Z ∗ ) exist under the following conditions:

X ∗ < K, ω3 > cD.

(iii) Stability properties of Model 3 ˜ EXY (X, ˆ Yˆ , 0) Feasibility and stability conditions of E0 (0, 0, 0), EX (K, 0, 0), EZ (0, 0, Z), ¯ and EXZ (K, 0, Z) of system (4) are given in Table 5. The coexistence equilibrium of Model (4) is E∗ (X ∗ , Y ∗ , Z ∗ ), where h i ω1 a1 X ∗ ∗ ∗ (Y + γ) − a2 + KωY ∗ (K − X ) ω3 − cD h i , Z∗ = Y∗ = ∗ c ω − β − a + ω1 a1 X (K − X ∗ ) 2

2

KωY ∗

such that K > X ∗ and ω3 > cD and X ∗ is the positive root of the cubic equation f (X ∗ ) = X ∗ 3 + (j − K)X ∗ 2 + j(a − K)X ∗ +

jωK(ω3 − cD) − jaK = 0. ca1

3 −cD) We note that 0 < X ∗ < K. Also, we have f (0) < 0 if ωKj(ω < jaK and f (K) = ca1 jωK(ω3 −cD) > 0 as ω3 > cD. Since f (0)f (K) < 0, there is a positive root of this equation lies ca1 3 −cD) in (0, K) when jωK(ω < jaK. Therefore, the non-zero equilibrium point E∗ (X ∗ , Y ∗ , Z ∗ ) ca1 exist under the following conditions:

X ∗ < K, ω3 > cD,

ωKj(ω3 − cD) < jaK. ca1

10

(iv) Stability properties of Model 4 ˜ EXY (X, ˆ Yˆ , 0) Feasibility and stability conditions of E0 (0, 0, 0), EX (K, 0, 0), EZ (0, 0, Z), ¯ and EXZ (K, 0, Z) of system (5) are given in Table 6. The coexistence equilibrium point of Model (5) is E∗ (X ∗ , Y ∗ , Z ∗ ), where     2 X∗ ∗ ∗ ∗ (δY + γ) a2 j + X + a − ω1 X ω3 − cD ∗ ∗      Y = , Z = c X ∗2 X ∗2 ∗ ∗ ∗ ∗ (ρY + β) ω1 X − a2 j + X + a − ω2 j + X + a such that K > X ∗ and ω3 > cD and X ∗ is the positive root of the cubic equation f (X ∗ ) = X ∗ 3 + (j − K)X ∗ 2 + j(a − K)X ∗ +

jωK(ω3 − cD) − jaK = 0 ca1

3 −cD) We note that 0 < X ∗ < K. Also, we have f (0) < 0 if ωKj(ω < jaK and f (K) = ca1 jωK(ω3 −cD) > 0 as ω3 > cD. Since f (0)f (K) < 0, there is a positive root of this equation lies ca1 3 −cD) in (0, K) when jωK(ω < jaK. Therefore, the non-zero equilibrium point E∗ (X ∗ , Y ∗ , Z ∗ ) ca1 exist under the following conditions:

X ∗ < K, ω3 > cD,

ωKj(ω3 − cD) < jaK. ca1

(v) Stability properties of Model 5: ˜ EXY (X, ˆ Yˆ , 0) Feasibility and stability conditions of E0 (0, 0, 0), EX (K, 0, 0), EZ (0, 0, Z), ¯ of system (6) are given in Table 7. The coexistence equilibrium point of and EXZ (K, 0, Z) Model (6) is E∗ (X ∗ , Y ∗ , Z ∗ ), where Y∗ =

(δY ∗ + γ1 )[a2 (γ + βY ∗ ) + X ∗ (a2 − ω1 )] ω3 − cD , Z∗ = c ω2 (X ∗ + βY ∗ + γ) − (ρY ∗ + β1 )[X ∗ (a2 − ω1 ) + a2 (γ + βY ∗ )]

such that K > X ∗ and ω3 > cD and X ∗ is the positive root of the equation     β(ω3 − cD) (a1 β − ω)(ω3 − cD) ∗ ∗2 X + (γ − K) + X −K γ+ = 0. c a1 c We note that 0 < X ∗ < K. Also, we have f (0) < 0 if aω1 < β and f (K) = ωK(ωc3 −cD) > 0 as ω3 > cD. Since f (0)f (K) < 0, there is a positive root of this equation lies in (0, K) when ω < β. Therefore, the non-zero equilibrium point E∗ (X ∗ , Y ∗ , Z ∗ ) exist under the following a1 conditions: ω X ∗ < K, ω3 > cD, < β. a1 • Stability result of E∗ (X∗ , Y∗ , Z∗ ) : Now we calculate the matrix A|τ =0 and B|τ =0 at the coexistence equilibrium point E∗ (X ∗ , Y ∗ , Z ∗ ), hence, by Routh-Hurwitz criterion, the 11

coexistence equilibrium point E∗ (X ∗ , Y ∗ , Z ∗ ) of systems (2)-(6) is asymptotically stable as all the roots of equation (10) have negative real parts when P + S, Q + T, R + U > 0 and (P + S)(Q + T ) > R + U . Due to the different parametric expressions of P, Q, R, S, T and U , we also get different parametric domains for the stability of the coexistence equilibrium point E∗ for systems (2)-(6). To represent any explicit result for each models are difficult due to the highly nonlinear mathematical expressions. So, we show the local stability results of the coexistence equilibrium point E∗ (X ∗ , Y ∗ , Z ∗ ) for systems (2)-(6) numerically by phase-space diagram (Figure 1(a)= system (2), Figure 1(b)= system(3), Figure 1(c)= system(4), Figure 1(d)= system(5), Figure 1(e)=system(6).

3.2

Stability properties of the delayed systems (1-5)

Now we establish the stability properties of systems (2)-(6) in generalized form. For so, we consider the linearized system (8) at coexistence equilibrium point E∗ (X ∗ , Y ∗ , Z ∗ ). If P + S, Q + T, R + U > 0 and (P + S)(Q + T ) > R + U happen for τ = 0, then we come forward to the next step of analysis. For the purpose of study the Hopf bifurcation, we take eigenvalues of characteristic equation (9) of systems (2)-(6) at coexistence equilibrium point E∗ (X ∗ , Y ∗ , Z ∗ ) as λ = iω, (ω > 0). Putting this value in equation (9) and then separating the real and imaginary parts, we get P ω 2 − R = (U − Sω 2 ) cos ωτ + T ω sin ωτ, ω 3 − Qω = T ω cos ωτ − (U − Sω 2 ) sin ωτ.

(11)

Now squaring and adding the both equations of (11), we get ω 6 + ξ1 ω 4 + ξ2 ω 2 + ξ3 = 0,

(12)

where ξ1 = P 2 − S 2 − 2Q, ξ2 = Q2 − T 2 + 2SU − 2P R, ξ3 = R2 − U 2 . Equation (12) has at least one real positive root. Without loss of generality, assume that it has six positive real roots, defining by ωl , (l = 1, 2, 3, 4, 5, 6) respectively. Now from equation (11), we get     (U − Sω 2 )(P ω 2 − R) + T ω 2 (ω 2 − Q) 1 −1 k cos + 2kπ , l = 1, ..., 6, k = 0, 1, ... τl = ωl (U − Sω 2 )2 + T 2 ω 2

then ±iωl are a pair of purely imaginary roots of equation (9) with τlk . Let us define τ0 = τl00 = min{τl0 }, ω0 = ωl0 . We have the following transversality conditions   −1 dλ ρ1 σ1 + ρ2 σ2 Re = > 0, dτ ρ21 + ρ22 τ =τ0 if σ2 > 0 and ρ1 , σ1 are of same sign, i.e., either ρ1 > 0, σ1 > 0 or ρ1 < 0, σ1 < 0, where ρ1 σ1 ρ2 σ2

= = = =

(U − Sω02 ), (−3ω02 + Q) cos ω0 τ0 − 2P ω0 sin ω0 τ0 + T, T ω0 , (−3ω02 + Q) sin ω0 τ0 − 2P ω0 cos ω0 τ0 + 2Sω0 . 12

Then, by the general Hopf bifurcation theorem, we have the following result on the stability and bifurcation of the systems (2)-(6). Theorem 3.3 If P + S, Q + T, R + U > 0 and (P + S)(Q + T ) > R + U happen when τ = 0, then assume τ 6= 0, ω0 is the unique positive root of equation (12), then the equilibrium E∗ (X ∗ , Y ∗ , Z ∗ ) is locally asymptotically stable for τ < τ0 and unstable for τ > τ0 . Furthermore, the systems (2)-(6) undergoes a Hopf-bifurcation at E∗ (X ∗ , Y ∗ , Z ∗ ) when τ = τ0 provided ρ1 σ1 + ρ2 σ2 > 0.

4 4.1

Finite Time Blow-Up Motivation, Control Mechanisms and “Ecological” Damping

In this section, we present results on finite time blow-up of the five delayed systems (2)-(6) and their non-delayed counterparts. From an ecological viewpoint, it is critical to derive the sharp conditions of population explosion/blow-up for any ecological model system. Parshad et al. [34] established that the solution to the three species model studied by Aziz-Alaoui [35], can blow-up in finite time even under the restriction derived in [35]. It is also useful then to study ecological mechanisms that can prevent such explosive phenomenon [30]. Recently, blow-up behavior is studied by many researchers for different dynamical systems. Zhang and Yang [36] obtained the conditions under which the solutions may exist globally or blow up in a finite time for a nonlinear parabolic equation subject to mixed boundary condition. Ling and Wang [37] studied blow-up criteria and global boundedness of non-negative solutions using super/sub solution method techniques. The exact conditions of blow-up and global existence of a nonlinear wave equation is studied by Jiang and Zhang [38]. Also, the sufficient conditions for the blow-up and global solutions are presented for nonlinear parabolic equation with different kinds of boundary conditions studied by Zhang and He [39]. Zhou et al. [40] illustrated that the blow-up behaviors of differential equations with piecewise constant arguments are quite different from those of the corresponding ordinary differential equations. Definition 1. Given a ODE model of a nonlinear process say, du = f (u), dt

(13)

we say finite time blow-up occurs if, lim

t→T ∗ <∞

kuk → ∞,

(14)

where k · k is the standard sup norm on Rn , u is the state variable in question that depends on time (t) and T ∗ is the blow-up time. In the context of population biology, finite time blow-up has also been well investigated [41, 42, 43, 44]. Note, we have now introduced an alternate viewpoint: finite time blow-up, 13

can be viewed as mimicing the explosive growth of an invasive species. This is formalised by equating: finite time blow-up = uncontrollable and unmanageable population level.

(15)

Here, the blow-up time T ∗ is viewed as the disaster time for the ecosystem. This is seen in many ecosystems worldwide, where large scale damages have been caused by invasive species [45, 46, 47, 48]. The blow-up mechanism has also been used to look at finite time extinction problems in Biology. Such as if a certain pest is desired to be eradicated in a finite amount of time [49, 50, 51]. Remark 1. Although populations cannot reach infinite values in finite time, they can grow rapidly [27]. For example, experimental evidence suggest that the human population may be growing hyperbolically, rather than logistically [52]. Data on the Burmese python suggests, that its population is growing at least exponentially [25]. Our approach investigates biological control mechanisms, that attempt to lower and control the targeted population before time T ∗ . This approach has distinct advantages: (i) There is no ambiguity as to what is a disastrous high level of population. (ii) There is a clear demarcation between when or if the disaster occurs. (iii) Our controls focus on avoiding classical chemical and biological controls. (iv) This method provides a predictive modeling tool for various ecological settings.

4.2

Finite time blow-up in the non delayed model system

Since all of the presented models have a similar generic structure, we prove the blow-up results for one case. We pick (3). We will first show that non-delayed model blows-up in finite time. Consider X wX 2 Y dX = a1 X(1 − ) − 2 , dt K X +b dY w1 X 2 Y w2 Y Z = −a2 Y + 2 − , dt X +b ρY Z + δY + β1 Z + γ1 dZ w3 Z 2 = cZ 2 − , dt Y +D

(16)

with suitable positive initial condition (X0 , Y0 , Z0 ). Theorem 1. Consider the model system given by (16), for any choice of parameters, and a δ1 > 0, such that c > δ1 , There exists initial data such that if this data meets the largeness condition ! (a2 + wβ12 ) |Y0 | w3 |Z0 | ln w3 > , |Y0 | > − D, (17) −D δ1 c − δ1 c−δ1 14

then the state variable Z will blow-up in finite time, that is lim

t→T ∗ <∞

kZk → ∞.

(18)

1 δ1 |Y0 |

Here the blow-up time T ∗ ≤

Proof. Consider the equation for the top predator dZ dt

= cZ 2 −

w3 Z 2 , Y +D

blow-up is trivial if
c > wD3 [31]. However if c < wD3 , in order to guarantee blow-up, we must 3 have that c − Yw+D > δ1 or equivalently we must guarantee that w3 Y > − D. (19) c − δ1 To this end, we will work with the equation for middle predator dY dt

w1 X 2 Y X 2 +b

= −a2 Y +

−

w2 Y Z , ρY Z+δY +β1 Z+γ1

We have via positivity dY dt

≥ −a2 Y −

w2 Y Z ρY Z+δY +β1 Z+γ1

≥ −a2 Y −

w2 Y Z β1 Z

= −(a2 +

w2 )Y, β1

we now use this to yield the following lower estimate on Y |Y | > |Y0 |e

w

−(a2 + β 2 )t 1

.

Note, we want equation (19) to hold. Thus it is sufficient to require |Y | > |Y0 |e To this end we consider

w

−(a2 + β 2 )t 1

w

−(a2 + β 2 )t

|Y0 |e

1

>

>

w3 c−δ1

w3 c−δ1

w (a2 + β 2 )t 1

− D.

− D.

We multiply both sides by e , and divide both sides by sides, and finally divide both sides by (a2 + wβ12 ) to yield,   |Y0 | 1 ln w3 −D > t. w (a2 + 2 ) β1

Herein we require

|Y0 | w3 −D c−δ1

− D, then take ln of both

c−δ1

> 1 for the ln term to stay positive. Note that dZ dt

blows-up at time T ∗ =

w3 c−δ1

= δ1 Z 2 ,

1 , δ1 |Z0 |

thus if we choose data such that   |Y0 | 1 1 , |Y0 | > ln w3 −D > t > T ∗ = δ1 |Z w (a2 + 2 ) 0| β1

c−δ1

w3 c−δ1

− D.

w3 Then the above guarantees that Y will remain above the critical level c−δ −D, for sufficiently 1 long enough time, for Z to blow-up. This yields that as long as the following holds   w (a2 + 2 ) |Y0 | w3 |Z0 | ln w3 −D > δ1β1 , |Y0 | > c−δ −D 1 c−δ1

15

For the results of the rest models, we follow the same process and the results are listed in Table(8) (non-delayed and delayed both of the cases).

4.3

Finite time blow-up in the delayed model system

Now, we consider the delayed model system (3): dX dt

= a1 X(1 −

dY dt

= −a2 Y +

dZ dt

= cZ 2 −

X ) K

−

wX 2 Y X 2 +b

,

w1 (X(t−τ ))2 Y (t−τ ) (X(t−τ ))2 +b

−

w2 Y Z , ρY Z+δY +β1 Z+γ1

w3 Z 2 . Y +D

Next, we show finite time blow-up in the delayed model systems via the following theorem. Theorem 2. Consider the delayed model system given by system (3), for any choice of parameters, and a δ1 > 0, such that c > δ1 , There exists initial data such that if this data meets the largeness condition ! (a2 + wβ12 ) |Y0 | , (20) |Z0 | ln w3 > δ1 −D c−δ1 then the state variable Z will blow-up in finite time, that is limt→T ∗ <∞ kZk → ∞. The proof is an easy extension of the non delayed case. Proof. Essentially, we have via positivity [33] of solutions to delayed system, dY dt

= −a2 Y +

w1 (X(t−τ ))2 Y (t−τ ) (X(t−τ ))2 +b

w2 Y Z ρY Z+δY +β1 Z+γ1

−

≥ −(a2 +

w2 )Y, β1

we now use this to yield the following lower estimate on Y w

−(a2 + β 2 )t

|Y | > |Y0 |e

1

and now proceed exactly as we did earlier to derive that the state variable Z will blow-up in finite time. Since the structure of the other model systems is the same, that is the delay terms appear in the middle predator equation, the same technique enables us to state the following corollary. Theorem 3. Consider the delayed model systems (2)-(6), for sufficiently large initial data the state variable Z for each model will blow-up in finite time, that is lim

t→T ∗ <∞

kZk → ∞. 16

Remark 2. Note, although our results are true for large initial data. There are several new results [53] that point to blow-up, when the initial data is small as well. It would make for extremely interesting work, if we could ascertain the effect of time delay on such systems. For example one might ask if time delay hinders or excerbates blow-up. Furthermore, it would also be extremely interesting to attempt to prove small data blow-up results, under a time delay. Remark 3. Note, with delay in the top predator, particularly in the sexual reproduction term, cZ 2 , global existence can be proved for any initial condition, and blow-up is prevented [54]. Thus the species in which gestation is the most relevant, seems critical to prevent population explosions.

5

Numerical simulation

To demonstrate the analytic results of models (2)-(6) and their non-delayed counter part, we choose five different parameter sets for each model which are given in Table 9. For each model, we show the stability of the coexistence equilibrium point (E∗ ) by phase-space diagram (Figure 1.n(i), n=a,b,c,d,e) and bifurcation diagram (Figure 1.n(ii), n=a,b,c,d,e). In Figure 1.n(i) (n=a,b,c,d,e), if we choose τ = 0 (black trajectory) and τ < τ0 (green trajectory), then the coexistence equilibrium point (E∗ ) is asymptotically stable and if we choose τ > τ0 , then the coexistence equilibrium point (E∗ ) is unstable and then solutions of the systems converge to the nearest stable limit cycle arround E∗ (blue trajectory). This dynamical variation occurs by changing the value of τ from 0 to higher value through its critical value τ0 . This dynamical phenomenon is also depicted by the bifurcation diagram of the systems (1-5) with the variation of τ in Figure 1.n(ii) (n=a,b,c,d,e). By Figure 1.n(ii) (n=a,b,c,d,e), we observe that switching of stability of the systems (1-5) is possible if we varry the parameter τ through its critical values and system is stable and unstable according as τ < τ0 and τ > τ0 .

5.1

Numerical simulation of blow-up

In this section, we explore the effect of the lag/delay on blow-up. We primarily ask the following questions 1. What is the effect of a time lag/delay on the dynamics of a system, that blows-up in finite time without delay? 2. Can we compare across models to see how various parameters influence blow up? 3. Based on the above, can we view delay as providing damping? In order to answer the first question, first we consider the system (2) with parameter set written in Table 10 and initial condition IC: (5, 6, 5), we see that system (2) without delay (figure 2(a)) and with delay τ = 0.231 (figure 2(b)) will go to steady state. Next we consider system (3) (parameter set is written in Table 10 and initial condition IC: (5, 6, 5)), the 17

w3 ) and figure D 3(b): with delay τ = 0.0231). System (4) does not blow-up in finite time for the IC: (5, 6, 5) (parameter set is in Table 10), without (figure 4(a)) or wit delay (figure 4(b)).However if we analogously view system (5) with parameter set as in Table 10 and IC: (5, 6, 5), note w3 c< , then we see the system will blow-up in finite time (see figure 5(a)) without delay. D Now we assume there is a delay, choosing τ = 1.9231, starting from the same IC: (5, 6, 5), we do not have blow-up up (see figure 5(b). The exact situation is seen with the system (6) in figure 6 respectively. Note, our goal here is to compare across models, the effects of delay and parameters on the phenomenon of finite time blow-up. To this end we start the delayed and non-delayed models from the same IC, in the simulations. system will blow-up in finite time (see figure 3(a): without delay (Note c <

In Figure (7), we observe the dynamical changes of each system (2)-(6) with respect to the varying initial conditions Y0 and Z0 along parameter set as in Table (10), where stable and blow up regions are stated as AS (Asymptotic Stable) and BL (Blow up) respectively. Black solid line in the Y0 Z0 -plane separates these two regions. Remark 4. Our results indicate that delay can certainly act as a damping mechanism if the gestation period, is in some critical range. However this certainly depends on the parameter ranges in question, as well as initial conditions. Our results enable us to answer the third question in the affirmative too. That is, delay can certainly provide damping, by preventing the top predator populations from exploding. It is interesting to explore this further in the context of various realistic gestation periods, and specific species.

6

Comparative study between different models

From the model formulation it is clear that the common issues of all the models are prey (X) grows logistically and top predator (Z) grows by sexual reproduction (for so mating of male and female happen and represented by Z × Z = Z 2 term in third equation), also Z is generalist type and it is represented by modified Leslie-Grower scheme (Third equation of each model). Only differences between the systems (2)-(6) come out with respect to the different combinations of functional responses of intermediate predator (Y ) and top predator (Z). It is also significant constituents of the systems (2)-(6) of that the functional response of Z is always complex than the functional response of Y which we mentioned in introduction. All the models (2)-(6) have five boundary equilibrium points and one coexistence equilibrium point. Feasibility and stability analysis of five boundary equilibrium points (E0 (0, 0, 0), ˜ EXY (X, ˆ Yˆ , 0) and EXZ (K, 0, Z)) ¯ of all the systems (2)-(6) in Section EX (K, 0, 0), EZ (0, 0, Z), 3, we see that these equilibrium points are all neutral point from the dynamical point of view (one eigenvalue of the Jacobian matrix at each boundary equilibrium point is zero). E0 (0, 0, 0) and Ex (K, 0, 0) are always feasible for all the systems. But feasibility conditions of ˜ Exy (ˆ x, yˆ, 0) for system (2)-(3) and system (4)-(5) remain same. Also feasibility of EZ (0, 0, Z) ¯ and EXZ (K, 0, Z) are same for all the systems. From the mathematical analysis, it is very clear and easy to say that these feasibility results depends upon the choices of the functional 18

responses of first predator upon the prey population. Also, it is very interesting to see from the analysis of coexistence equilibrium point E∗ (X ∗ , Y ∗ , Z ∗ ) that analytic value of Y ∗ for all the models are same which is coming from the third equations of all models. But expressions of X ∗ and Z ∗ of all the systems are different due to different combinations of functional responses for intermediate and top predator. For so, feasibility conditions emphasize different parametric domains which are significantly followed by the stability characterization of this equilibrium point. Trivial condition of finite time blow-up is c > wD3 i.e., same which is coming from the third equation which of all models. But the non-trivial conditions of finite time blow-up for all the models (non-delayed and as well as delayed) are different those are coming from the second equation of each model and significantly depending upon the functional response of top predator Z.

7

Conclusions and Discussions

In this paper, a general framework has been structured to study the dynamical behavior of hybrid tritrophic food chain model consisting specialist middle predator and sexually reproductive generalist top predator. The different combinations of prey dependent and preypredator dependent functional responses are used to formulate five differnt model systems. We discussed the stability properties of each equilibrium points for non-delayed couter-part of each model system. Existence of Hopf-bifurcation for coexisting equilibrium point of each delayed model system has been investigated under a general result. Blow-up phenomena in finite time has been observed for some of the model according to the presence of function response. Numerical simulations are presented to verify the analytical predictions. Note that the generalist predators are considered poor agents for biological control due to a tendency for frequent interference. However real data contradicts this [55], rather ”paradoxically”. Earlier works have investigated this very phenomenon [32]. What we now see is that interference of prey handling predators, is a critical factor in facilitating population explosions. That is if we compare models (3) to (2) and (5) to (4), we see that removing the interference of prey handling predators (that is setting ρ = 0) immediately damps the explosive population. However, it now seems that the feedback loop caused by interference might depend in an intricate way on time delay or gestation. This needs to be investigated in much more detail, in light of an appropriate damping [30] (here incorporated via gestation effect). To some degree this is accomplished in the current manuscript. We find from a numerical point of view is the evidence that gestation delay, can act as a damping mechanism to prevent population explosions. What is observed in our numerical simulations, is that for the set of parameters considered, there is a critical lower bound for the gestation delay in the middle predator, such that values larger than this act to prevent blow-up, but smaller values don’t. Ecologically, this says that there is possibly a complex feedback process between the gestation periods of the top and middle predator, that couple in a way to prevent explosion. It seems in the current model the gestation of middle predator 19

(which is prey for the top predator) is certainly more critical. It would be interesting to think of this from a biological control viewpoint, in particular in the advent of climate change. For example it has been noted that very recently the python population in the glades has declined (after increasing exponentially for several years). This has been attributed to unusually cold recent winters [56]. It would be interesting to investigate if such climate change can have an effect on the gestation periods of certain prey of the python, in particular the prey noted in [25]. The current results could then serve to explain this situation, if the gestation period increase, even slightly under such extreme environmental/climatic fluctuation. Since, these food chain models represent a wide range of real-world systems, the findings are interesting from the point of view of applications. In recent world, different aspects of nature are coming out as research objective. So researchers are exploring these theoretically and also validating with real world data. Real world applications are paying a strong insight in science and technology. Application of stochastic process on delayed neural network model [57], nonfragile asynchronous control for uncertain chaotic Lurie network systems with Bernoulli stochastic process [58], stability analysis of fuzzy sampled-data Markovian chaotic systems [59] are genuine examples for applicable mathematics with contemporary outlook. Similarly in our nature lots of examples are there which is shown finite time blow up. Several biological conditions will be coming out to explore the biological phenomenon of invasive occurrence which will be followed by new techniques of Mathematics.

Acknowledgments The authors are grateful to the anonymous referees and Editor in Chief of this journal for their careful reading, valuable comments and helpful suggestions, which have helped them to improve the presentation of this work significantly.

References [1] Wang, X., Peng, M., Liu, X.: Stability and Hopf bifurcation analysis of a ratio-dependent predatorprey model with two time delays and Holling type III functional response. Applied Mathematics and Computation, 268 (2015), 496-508. [2] Sahoo, B., Poria, S.: Diseased prey predator model with general Holling type Interactions. Applied Mathematics and Computation, 226 (2014), 83-100. [3] Zhang, S., Wang, F., Chen, L.: A food chain model with impulsive perturbations and Holling IV functional response. Chaos, Solitons & Fractals, 26 (2005), 855-866. [4] Yin, H., Zhou, J., Xiao, X., Wen, X.: Analysis of a diffusive Leslie-Gower predator-prey model with nonmonotonic functional response. Chaos, Solitons & Fractals, 65 (2014), 51-61. 20

[5] Holling, C. S.: The functional response of predators to prey density and its role in mimicry and population dynamics. The Memoirs of the Entomological Society of Canada, 97(45)(1965), 560. [6] Beddington, J. R.: Mutual interference between parasites or predators and its effect on searching efficiency. Journal of Animal Ecology, 44(1)(1975), 331-340. [7] Crowley, P. H., Martin, E. K.: Functional responses and interference within and between year classes of a dragonfly population. Journal of the North American Benthological Society, 8(3)(1989), 211-221. [8] Hassell, M. P., Varley, G. C.: New inductive population model for insect parasites and its bearing on biological control. Nature, 223(1969), 1133-1137. [9] Cui, Q., Zhang, Q., Qiu, z., Hu, Z.: Complex dynamics of a discrete-time predator-prey system with Holling IV functional response. Chaos, Solitons & Fractals, 87(2016), 158171. [10] Arditi, R., Akcakaya, H. R.: Variation in plankton densities among lakes: a case for ratio-dependent predation models. The American Naturalist, 138(5)(1991), 1287-1296. [11] Gutierrez, A. P.: Physiological basis of ratio-dependent predator-prey theory: the metabolic pool model as a paradigm. Ecology, 73(5)(1992), 1552-1563. [12] Sarwardi, S., Haque, M., Mandal, P. K.: Persistence and global stability of bazykin predator-prey model with beddington-deangelis response function. Communication in Nonlinear Science and Numerical Simulation, 19(1)(2014), 189-209. [13] Tripathi, J. P., Abbas, S., Thakur, M.: Dynamical analysis of a prey-predator model with beddington-deangelis type function response incorporating a prey refuge. Nonlinear Dynamics, 80(2015), 1-20. [14] Skalski, G. T., Gilliam, J. F.: Functional responses with predator interference: viable alternatives to the holling type II model. Ecology, 82(11)(2011), 3083-3092. [15] Tucker, M. A., Rogers, T. L.: Examining predatorprey body size, trophic level and body mass across marine and terrestrial mammals. Proceedings of the Royal Society B-Biological Sciences, 281(1797) (2014), 20142103. [16] Agrawal, R., Jana, D., Upadhyay, R. K., Rao, V. S. H.: Dynamic relationship between mutual interference and gestation delay of hybrid tritrophic food chain model. ANZIAM, 59(3) (2018), 370-401. [17] Upadhyay, R. K., Chattopadhyay, J.: Chaos to order: role of toxin producing phytoplankton in aquatic systems. Nonlinear Anal Modell Control, 10(4) (2005), 383396. [18] Upadhyay, R. K., Iyengar, S. R. K., Rai, V.: Chaos: an ecological reality? Int J Bifurc Chaos, 8 (1998), 13251333.

21

[19] Upadhyay, R. K., Rai, V., Iyengar, S. R. K.: Species extinction problem: genetic vsecological factors. Appl Math Modell, 25 (2001), 937951. [20] Jana, D., Tripathi, J. P.: Impact of generalist type sexually reproductive top predator interference on the dynamics of a food chain model. International Journal of Dynamics and Control, 5 (2017), 9991009. [21] Parshad, R. D., Basheer, A., Jana, D., Tripathi, J.: Do prey handling predators really matter: Subtle effects of a Crowley-Martin functional response. Chaos, Solitons and Fractals, 103 (2017), 410-421. [22] Jana, D., Agrawal, R., Upadhyay, R. K.: Dynamics of generalist predator in a stochastic environment: effect of delayed growth and prey refuge. Applied Mathematics and Computation, 268 (2015), 1072-1094. [23] Jana, D., Gopal, R., Lakshmanan, M.: Complex dynamics generated by negative and positive feedback delays of a preypredator system with prey refuge: Hopf-bifurcation to Chaos. International Journal of Dynamics and Control, 5(4) (2017), 1020-1034. [24] Parshad, R. D., Kumari, N., Kouachi, S.: A remark on ”Study of a Leslie-Gower-type tritrophic population model. Chaos, Solitons and Fractals, 14 (2002), 1275-1293. Chaos, Solitons & Fractals, 71(2015) 22-28. [25] Dorcas, M. E., Willson, J. D., Reed, R. N., Snow, R. W., Rochford, M. R., Miller, M. A., Mehsaka, W. E., Andreadis, J. P. T., Mazzotti, F. J., Romagosa, C. M., Hart, K. M.: Severe Mammal Declines Coincide with Proliferation of Invasive Burmese Pythons in Everglades National Park. Proceedings of the National Academy of Sciences, 109 (2012), 2418-2422. [26] Pimentel, D., Zuniga, R., Morrison, D.: Update on the environmental and economic costs associated with alien-invasive species in the united states. Ecological Economics, 52(3) (2005), 273-288. [27] Berryman, A.: The theory and classification of outbreaks. Insect Outbreaks, Academic Press, San Diego, CA, 1987. [28] Ludwig, D., Jones, D., Holling, C.: Qualitative analysis of insect outbreak systems: The spruce budworm and forest. The journal of animal ecology, 47(1)(1978), 315-332. [29] Driesche, R. V., Bellows, T.: Biological Control. Kluwer Academic Publishers, Massachusetts, 1996. [30] Parshad, R. D., Qansah, E., Black, K., Beauregard, M.: Biological control via ”ecological” damping: An approach that attenuates non-target effects. Mathematical Biosciences, 273 (2016), 23-44. [31] Parshad, R. D., Abderrahmanne, H., Upadhyay, R. K., Kumari, N.: Finite time blowup in a realistic food chain model, ISRN Biomathematics, 2013 (2013), Article ID 424062. 22

[32] Parshad, R. D., Bhowmick, S., Quansah, E., Basheer, A., Upadhyay, R. K.: Predator interference effects on biological control: the paradox of the generalist predator revisited. Communications in Nonlinear Science and Numerical Simulation, 39 (2016), 169-184. [33] Hale, J.: Theory of Functional Differential Equations. Springer-Verlag, Berlin, 1977. [34] Parshad, R. D., Kumari, N., Kouachi, S.: A remark on ”Study of a Leslie-Gower-type tritrophic population model”. Chaos, Solitons & Fractals, 71 (2015), 22-28. [35] Aziz-Alaoui, M. A.: Study of a Leslie-Gower type tri-trophic population model. Chaos Solitons & Fractals, 14(8) (2002), 1275-1293. [36] Zhang, R., Yang, Z.: Uniform blow-up rates and asymptotic estimates of solutions for diffusion systems with weighted localized sources. Journal of Applied Mathematics and Computation, 32 (2010), 429-441. [37] Ling, Z., Wang, Z.: Global boundedness and blow-up for a parabolic system with positive Dirichlet boundary value. Journal of Applied Mathematics and Computation, 46 (2014), 123-134. [38] Jiang, Y., Zhang, Y.: Exact conditions of blow-up and global existence for the nonlinear wave equation with damping and source terms. Nonlinear Dynamics, 76(1) (2014), 139146. [39] Zhang, T., He, Y.: Blow-up and global solutions for a class of nonlinear parabolic equations with different kinds of boundary conditions. Applied Mathematics and Computation, 217 (2010), 801-810. [40] Zhou, Y. C., Yang, Z. W., Zhang, H. Y., Wang, Y.: Theoretical analysis for blow-up behaviors of differential equations with piecewise constant arguments. Applied Mathematics and Computation, 274 (2016), 353-361. [41] Kim, K., Lin, Z.: Blow-up in a three species cooperating model. Applied Mathematics Letters, 17 (2004), 89-94. [42] Lou, Y., Nagylaki, T., Ni, W.: On diffusion induced blowups in a mutualistic model. Nonlinear Analysis, 45 (2001), 329-342. [43] Lou, Y., Munther, D.: Dynamics of a three species competition model. Discrete & Continuous Dynamical Systems-A, 32 (2012), 3099-3131. [44] Hillen, T., Painter, K.: A users guide to PDE models for chemotaxis. Journal of Mathematical Biology, 57 (2009), 183-217. [45] Seebens, H., Blackburn, T. M., Dyer, E. E., Genovesi, P., Hulme, P. E., Jeschke, J. M., Bacher, S.: No saturation in the accumulation of alien species worldwide. Nature Communications, 8 (2017), 14435. [46] Lewis, M. A., Petrovskii, S. V., Potts, J. R.: The mathematics behind biological invasions. Springer, 44 (2016). 23

[47] Dorcas, M. E., Willson, J. D., Reed, R. N., Snow, R. W., Rochford, M. R., Miller, M. A., Mehsaka, W. E., Andreadis, P. T., Mazzotti, F. J., Romagosa, C. M., Hart, K. M.: Severe Mammal Declines Coincide with Proliferation of Invasive Burmese Pythons in Everglades National Park. Proceedings of the National Academy of Sciences, 109 (2012), 2418-2422. [48] Letnic, M., Webb, J., Shine, R.: Invasive cane toads (Bufo marinus) cause mass mortality of freshwater crocodiles (Crocodylus johnstoni) in tropical Australia. Biological Conservation, 141 (2008), 17731782. [49] Parshad, R. D., Wickramsooriya, S., Bailey, S.: A remark on ”Biological control through provision of additional food to predators: A theoretical study” [Theoretical Population Biology 72 (2007) 111-120]. In Revision, Theoretical Population Biology (2019). [50] Lundgren, J. G., Fergen, J. K.: Predator community structure and trophic linkage strength to a focal prey. Molecular ecology, 23(15) (2014), 3790-3798. [51] Sappington, T.: Emerging issues in Integrated Pest Management implementation and adoption in the North Central USA. In Integrated Pest Management (pp. 65-97). Springer, Dordrecht (2014). [52] Grinn, L., Hermann, P., Korotayev, A., Tausch, A.: History & Mathematics: Processes and Models of Global Dynamics. Volgograd ’Uchitel’ Publishing House, 2010. [53] Parshad, R. D., Quansah, E., Beauregard, M., Kouachi, S.: On ”small” data blow-up in a three species food chain model. Computers and Mathematics with Applications, 73(2) (2017), 576-587. [54] Parshad, R.D., Upadhyay, R.K., Mishra, S., Tiwari, S., Sharma, S.: On the explosive instability in a three species food chain model with modified Holling type IV functional response. Mathematical Methods in the Applied Sciences, (2017), DOI: 10.1002/mma.4419, Appeared online May 3rd 2017. [55] Rosenheim, J. A., Kaya, H. K., Ehler, L. E., Marois, J. J., Jaffee, B. A.: Intraguild predation among biological-control agents: Theory and practice. Biological Control, 5 (1995), 303-335. [56] Mazzotti, F. J., Cherkiss, M. S., Hart, K. M., Snow, R. W., Rochford, M. R., Dorcas, M. E., Reed, R. N.: Cold-induced mortality of invasive burmese pythons in south orida. Biological Invasions, 13(1) (2011), 143-151. [57] Wang, J., Shi, K., Huang, Q., Zhong, S., Zhang, D.: Stochastic switched sampled-data control for synchronization of delayed chaotic neural networks with packet dropout. Applied Mathematics and Computation, 335 (2018), 211-230. [58] Shi, K., Tang, Y., Zhong, S., Yin, C., Huang, X., Wang, W.: Nonfragile asynchronous control for uncertain chaotic Lurie network systems with Bernoulli stochastic process. International Journal of Robust and Nonlinear Control, 28 (2018), 1693-1714. 24

[59] Zhang, R., Liu, X., Zeng, D., Zhonga, S., Shi, K.: A novel approach to stability and stabilization of fuzzy sampled-data Markovian chaotic systems. Fuzzy Sets and Systems, 344 (2018), 108-128.

25

Table 1: Functional form of system (2)-(6) Model

f1 (X, Y, 0)

(2)

a1 X(1 −

X ) K

−

wX 2 Y X 2 +b

(3)

a1 X(1 −

X ) K

−

wX 2 Y X 2 +b

(4)

a1 X(1 −

X ) K

−

X2 j

(5)

a1 X(1 −

X ) K

−

X2 j

(6)

a1 X(1 −

X ) K

−

f¯2 (X(t − τ ), Y (t − τ ), 0)

f2 (0, Y, Z)

w2 Y Z Y +βZ+γ

w1 (X(t−τ ))2 Y (t−τ ) (X(t−τ ))2 +b

cZ 2 −

w3 Z 2 Y +D

w2 Y Z ρY Z+δY +β1 Z+γ1

w1 (X(t−τ ))2 Y (t−τ ) (X(t−τ ))2 +b

cZ 2 −

w3 Z 2 Y +D

w2 Y Z Y +βZ+γ

w1 X(t−τ )Y (t−τ ) X 2 (t−τ ) +X(t−τ )+a j

cZ 2 −

w3 Z 2 Y +D

w2 Y Z ρY Z+δY +β1 Z+γ1

w1 X(t−τ )Y (t−τ ) X 2 (t−τ ) +X(t−τ )+a j

cZ 2 −

w3 Z 2 Y +D

w2 Y Z ρY Z++δY +β1 Z+γ1

w2 X(t−τ )Y (t−τ ) X(t−τ )+βY (t−τ )+γ

cZ 2 −

w3 Z 2 Y +D

−a2 Y −

−a2 Y −

wXY +X+a

wXY +X+a

w2 XY X+βY +γ

−a2 Y −

−a2 Y −

−a2 Y −

Table 2: Entities of Jacobian matrix for system (2)-(6) Model

Jacobian Matrix entities 2X wX 2 2wbXY ) − (X 2 +b)2 , f1Y = − X 2 +b , K (βZ+γ)w2 Z (Y +γ)w2 Y f2Y = −a2 − (Y +βZ+γ)2 , f2Z = − (Y +βZ+γ)2 , 2w1 bXY ¯ w1 X 2 f¯2X = (X 2 +b)2 , f2Y = X 2 +b , 2 w3 3Z f3Y = (Yw+D) 2 , f3Z = 2Z(c − Y +D ). 2wbXY wX 2 f1X = a1 (1 − 2X ) − (X 2 +b)2 , f1Y = − X 2 +b , K (β1 Z+γ1 )w2 Z (δY +γ1 )w2 Y −a2 − (ρY Z+δY , f2Z = − (ρY Z+δY +β 2, +β1 Z+γ1 )2 1 Z+γ1 ) 2 2w1 bXY ¯ w1 X ¯ f2X = (X 2 +b)2 , f2Y = X 2 +b , 2 w3 3Z f3Y = (Yw+D) 2 , f3Z = 2Z(c − Y +D ).

f1X = a1 (1 −

(2)

(3) f2Y =

2

(4)

f1X = a1 (1 −

2X ) K

f2Y = −a2 −

f1X = a1 (1 − f2Y = −a2 −

+X+a)2 (βZ+γ)w2 Z , f2Z (Y +βZ+γ)2

, f1Y = − X 2 wX =

2

2

, f¯2Y =

w3 Z 2 , f3Z (Y +D)2

2X ) K

−

2

( Xj +X+a)2 2

+X+a

,

w3 ). Y +D

, f1Y = − X 2 wX

, f2Z =

2

w1 Y (a− Xj )

w1 X

X2 j

,

2

+X+a)2

(ρY Z+δY +β1 Z+γ1 )2

+X+a j (Y +γ)w2 Y − (Y +βZ+γ) 2,

= 2Z(c −

wY (a− Xj )

2 ( Xj (β1 Z+γ1 )w2 Z

f¯2X =

, f¯2Y =

,

+X+a j (δY +γ1 )w2 Y − (ρY Z+δY +β 2, 1 Z+γ1 ) w1 X

X2 j

+X+a

,

w3 Z 3 , f3Z = 2Z(c − Yw+D ). (Y +D)2 (βY +γ)w Y (X+γ)w2 X 2 f1X = a1 (1 − 2X ) − , f = − , 1Y K X+βY +γ (X+βY +γ)2 (β1 Z+γ1 )w2 Z (δY +γ1 )w2 Y f2Y = −a2 − (ρY Z+δY +β Z+γ )2 , f2Z = − (ρY Z+δY +β Z+γ )2 , 1 1 1 1 (βY +γ)w2 Y ¯ (X+γ)w2 X f¯2X = (X+βY +γ) 2 , f2Y = X+βY +γ , 2 w3 3Z f3Y = (Yw+D) 2 , f3Z = 2Z(c − Y +D ).

f3Y =

(6)

2 ( Xj

( Xj +X+a)2

f3Y = (5)

wY (a− Xj )

−

w1 Y (a− Xj )

f¯2X =

f3 (0, Y, Z)

26

Table 3: Stability analysis of boundary equilibria of model system (2) Equilibrium and coordinate

Feasibility condition

E0 (0, 0, 0)

always

(i)

(ii)

EX (K, 0, 0)

˜ EZ (0, 0, Z)

(iii)

  

always

c=



ˆ Yˆ , 0) EXY (X,

Yˆ = (v)

q

a2 b a2 −w1 ˆ X ˆ 2 +b) a1 (K−X)( ˆ KwX

ˆ = X

¯ EXZ (K, 0, Z) ¯ is any +ve constant Z

−a1 0 0

2

wK −K 2 +b

2

Stability status Neutral 

 0  0

w1 K 2 , λ3 K 2 +b

Neutral

=0

0



˜ ˜  (β Z+γ)w 2Z  0 2 ˜ (β Z+γ)  2 ˜ w3 w3 Z ˜ 2Z(c − D ) D2 ˜ ˜ (β Z+γ)w w3 2Z ˜ λ1 = a1 , λ2 = −a2 − 2 , λ3 = 2Z(c − D ) ˜ (β Z+γ)  2 ˆ ˆ ˆ ˆ X X XY a1 (1 − 2K −w 0 ) − 2wb ˆ 2 +b)2 ˆ 2 +b (X X  ˆ +γ)w2 Y ˆ ˆY ˆ  (Y 2w1 bX 0 − ˆ  ˆ 2 +b)2 (X (Y +γ)2

a2 > w1 ,

−a2 −

0

ˆ K>X

c=

0

1K −a2 + w K 2 +b 0

λ1 = −a1 , λ2 = −a2 + a1 0

  0  0

w3 D

˜ is any +ve constant Z (iv)

Jacobian matrix  and eigenvalues  a1 0 0  0 −a2 0  0 0 0 λ1 = a1 , λ2 = −a2 , λ3 = 0

λ1 = a1 (1 −



w3 D

λ1 =

0

ˆ 2X ) K

−

ˆY ˆ 2wbX ˆ 2 +b)2 , λ2 (X

0

= 0, λ3 = 0

Neutral =0 

Neutral



Neutral

  

wK 2 −K 0 2 +b  ¯ ¯ (β Z+γ)w2 Z w1 K 2  + K 2 +b 0 −a2 − (β Z+γ) 0 2 ¯  ˜2 w3 Z w 3 ¯ 0 ) 2 Z(c − D D2 2 ¯ ¯ (β Z+γ)w 2Z 1K ¯ − w3 ) −a1 , λ2 = −a2 − (β Z+γ) +w , λ3 = 2Z(c 2 ¯ D K 2 +b

  

−a1

27

=0

Table 4: Stability analysis of boundary equilibria of model system (3) Equilibrium and coordinate

Feasibility condition

E0 (0, 0, 0)

always

(i)

(ii)

EX (K, 0, 0)

˜ EZ (0, 0, Z)

(iii)

ˆ Yˆ , 0) EXY (X,

Yˆ = (v)

q

a2 b a2 −w1 ˆ X ˆ 2 +b) a1 (K−X)( ˆ KwX

ˆ = X

¯ EXZ (K, 0, Z) ¯ is any +ve constant Z

 

always

c=

w3 D

˜ is any +ve constant Z (iv)



a2 > w1 ,

Jacobian matrix  and eigenvalues  a1 0 0  0 −a2 0  0 0 0 λ1 = a1 , λ2 = −a2 , λ3 = 0



−a1 0 0

2

wK −K 2 +b

c=

w3 D

Neutral 

 0  0

w1 K 2 , λ3 K 2 +b

Neutral

=0

0



˜ ˜  (β1 Z+γ 1 )w2 Z  0 2 ˜ (β1 Z+γ  1) 2 ˜ w3 Z w3 ˜ 2Z(c − D ) D2 ˜ ˜ )w Z (β Z+γ ˜ − w3 ) λ1 = a1 , λ2 = −a2 − 1 ˜ 1 22 , λ3 = 2Z(c D (β1 Z+γ1 )  2 ˆ ˆ ˆ ˆ X X XY a1 (1 − 2K −w 0 ) − 2wb 2 2 2 ˆ ˆ (X +b) X +b  ˆ +γ1 )w2 Y ˆ ˆY ˆ  (δ Y 2w1 bX 0 −  ˆ 2 +b)2 ˆ +γ1 )2 (X (δ Y

−a2 −

0

ˆ K>X

2

1K −a2 + w K 2 +b 0

λ1 = −a1 , λ2 = −a2 + a1 0

  0  0

0

Stability status

λ1 = a1 (1 −

   

0

ˆ 2X ) K

−

ˆY ˆ 2wbX ˆ 2 +b)2 , λ2 (X

2

wK −K 2 +b

−a1 0 0

λ1 = −a1 , λ2 =

0

= 0, λ3 = 0

0

Neutral =0 

Neutral



Neutral

  

 2 ¯ ¯ (β1 Z+γ 1 )w2 Z w1 K  −a2 − (β 0 2 ¯  K 2 +b 1 Z+γ1 ) ¯2 w3 Z w 3 ¯ ) 2 Z(c − D D2 2 ¯ ¯ (β1 Z+γ 1 )w2 Z 1K ¯ − w3 ) −a2 − (β +w , λ3 = 2Z(c 2 ¯ D K 2 +b 1 Z+γ1 )

28

=0

Table 5: Stability analysis of boundary equilibria of model system (4) Equilibrium and coordinate

Feasibility condition

(i)

E0 (0, 0, 0)

always

(ii)

EX (K, 0, 0)

always

Jacobian matrix  and eigenvalues  a1 0 0  0 −a2 0  0 0 0 λ1 = a1 , λ2 = −a2 , λ3 = 0   −a1 − K 2 wK 0 +K+a j    0 0  −a2 + K 2 wK   +K+a

Stability status Neutral

Neutral

j

˜ EZ (0, 0, Z)

(iii)

c=

w3 D

˜ is any +ve constant Z

ˆ Yˆ , 0) EXY (X,

(iv)

ˆ = X Yˆ =

∆ = (w1 − (v)

Kw a2 )2

ˆ K > X,

ˆ +X+a)

−

4aa2 2

(w1 − a2 )2 >

j

¯ EXZ (K, 0, Z)

¯ is any +ve constant Z

c=

w3 D

(β Z+γ)

 a1 (1 −     

w1 > a2 ,

√ (w1 −a2 )+ ∆

2a2 j ˆ2 ˆ X a1 (K−X)( j



0 0 0 λ1 = −a1 , λ2 = −a2 + K 2 wK , λ3 = 0 +K+a j   a1 0 0 ˜ ˜   (β Z+γ)w 2Z  0  −a2 − 0 2 ˜ (β Z+γ)   2 ˜ w3 Z w3 ˜ 0 2 Z(c − ) 2 D D ˜ ˜ (β Z+γ)w w3 2Z ˜ λ1 = a1 , λ2 = −a2 − 2 , λ3 = 2Z(c − D ) = 0 ˜ ˆ 2X ) K

ˆ2

−

ˆ (a− X ) wY j

ˆ2 ˆ +X+a

− Xˆ 2wX

ˆ2 2 ˆ ( Xj +X+a) ˆ2 X ˆ w1 Y (a− j ) ˆ2 2 ˆ ( Xj +X+a)

j

−

0

0

0 ˆ +γ)w2 Y ˆ (Y ˆ +γ)2 (Y

0

λ1 = a1 (1 −

ˆ 2X ) K

0

ˆ2

−

ˆ (a− X ) wY j ˆ2

x+a)2 ( Xj +ˆ

, λ2 = 0, λ3 = 0

Neutral

      

Neutral

4aa2 2 j

     

− K 2wK

−a1 0

−a2 −

0

λ1 = −a1 , λ2 = −a2 −

29

2

+K+a 2 ¯ ¯ (β Z+γ)w 2Z 1K + Kw 2 ¯ (β Z+γ)2 +K+a

0

j

¯2 w3 Z D2 ¯ ¯ (β Z+γ)w2 Z 2 ¯ (β Z+γ)

0

j

+

w1 K 2

K2 j

+K+a

¯ − w3 ) 2Z(c D ¯ − , λ3 = 2Z(c

     

w3 ) D

Neutral

=0

Table 6: Stability analysis of boundary equilibria of model system (5) Equilibrium and coordinate

Feasibility condition

(i)

E0 (0, 0, 0)

always

(ii)

EX (K, 0, 0)

always

Jacobian matrix  and eigenvalues  a1 0 0  0 −a2 0  0 0 0 λ1 = a1 , λ2 = −a2 , λ3 = 0   −a1 − K 2 wK 0 +K+a j    0 0  −a2 + K 2 wK   +K+a

Stability status Neutral

Neutral

j

˜ EZ (0, 0, Z)

(iii)

c=



(iv)

ˆ = X Yˆ =

2a2 j ˆ2 ˆ X a1 (K−X)( j

∆ = (w1 − (v)

Kw a2 )2

     

ˆ K > X,

ˆ +X+a)

−



w1 > a2 ,

√ (w1 −a2 )+ ∆

4aa2 2

(w1 − a2 )2 >

j

¯ EXZ (K, 0, Z)

¯ is any +ve constant Z

c=

w3 D

a1

0

K2 j

wK , λ3 +K+a

0

  0  0

w3 D

˜ is any +ve constant Z

ˆ Yˆ , 0) EXY (X,

0 0 λ1 = −a1 , λ2 = −a2 +

0

=0 

˜ ˜  (β1 Z+γ 1 )w2 Z  0 2 ˜ (β1 Z+γ  1) 2 ˜ w3 Z w3 ˜ 2 Z(c − ) 2 D D ˜ ˜ )w Z (β Z+γ ˜ − w3 ) λ1 = a1 , λ2 = −a2 − 1 ˜ 1 22 , λ3 = 2Z(c D (β1 Z+γ1 ) 2 ˆ ˆ (a− X ) wY ˆ2 ˆ j wX 2X − Xˆ 2 0 a1 (1 − K ) − Xˆ 2 2 ˆ ˆ +X+a) +X+a (

−a2 −

j

=0

j

ˆ2 ˆ (a− X w1 Y ) j

−

0

ˆ2 2 ˆ ( Xj +X+a)

0

0

λ1 = a1 (1 −

Neutral

ˆ 2X ) K

ˆ +γ1 )w2 Y ˆ (δ Y ˆ +γ1 )2 (δ Y

0

ˆ2

−

ˆ (a− X ) wY j ˆ2

2 ˆ ( Xj +X+a)

, λ2 = 0, λ3 = 0

      

Neutral

4aa2 2 j

     

− K 2wK

−a1 0

−a2 −

0

λ1 = −a1 , λ2 = −a2 −

30

2

+K+a 2 ¯ ¯ (β1 Z+γ 1 )w2 Z 1K + Kw 2 ¯ (β1 Z+γ1 )2 +K+a

0

j

¯2 w3 Z D2 ¯ ¯ (β1 Z+γ1 )w2 Z 2 ¯ (β1 Z+γ 1)

0

j

+

w1 K 2

K2 j

+K+a

¯ − w3 ) 2Z(c D ¯ − , λ3 = 2Z(c

     

w3 ) D

Neutral

=0

Table 7: Stability analysis of boundary equilibria of model system (6) Equilibrium and coordinate

Feasibility condition

(i)

E0 (0, 0, 0)

always

(ii)

EX (K, 0, 0)

always

(iii)

˜ EZ (0, 0, Z)

c=

w3 D

˜ is any +ve constant Z ˆ Yˆ , 0) EXY (X,

(iv)

ˆ = X

q −A2 + A2 2 −4A1 A3

Jacobian matrix  and eigenvalues  a1 0 0  0 −a2 0  0 0 0 λ1 = a1 , λ2 = −a2 , λ3 = 0   (K+γ)w2 K −a1 − 0 K+γ   (K+γ)w2 K  0 −a2 + 0  K+γ 0 0 0 (K+γ)w2 K λ1 = −a1 , λ2 = −a2 + , λ3 = 0 K+γ   a1 0 0 ˜ ˜   Z )w (β1 Z+γ 1 1   0 0 −a2 − 2 ˜ (β1 Z+γ   1) ˜2 w3 Z w 3 ˜ − ) 2Z(c 0 D D2 ˜ ˜ )w Z (β Z+γ ˜ − w3 ) = 0 λ1 = a1 , λ2 = −a2 − 1 ˜ 1 12 , λ3 = 2Z(c D

always,

,

   

a1 (1 −

0 λ1 = a1 (1 −

2A1 ˆ (w1 −a2 )X−a 2γ , β a1 w1 β A1 = wK > 0, 1 A2 = w1 − a2 − a1 βw , w

Yˆ =

A3 = −a2 γ < 0

(v)

¯ EXZ (K, 0, Z) ¯ is any +ve constant Z

c=



w3 D

λ1 =

(β1 Z+γ1 ) ˆ ˆ +γ)w2 Y ˆ ˆ ˆ (β Y (X+γ)w 2X 2X ) − − ˆ ˆ +γ ˆ ˆ +γ)2 K X+β Y (X+β Y ˆ +γ)w2 Y ˆ (β Y 0 ˆ ˆ +γ)2 (X+β Y

0 ˆ 2X ) K

−

ˆ +γ)w2 Y ˆ (β Y ˆ ˆ +γ)2 , λ2 (X+β Y

0 −

ˆ +γ1 )w2 Y ˆ (δ Y (δ y ˆ+γ1 )2

0

= 0, λ3 = 0

Stability status Neutral

Neutral

Neutral

   



(K+γ)w2 K − 0 K+γ  ¯ ¯ (β1 Z+γ1 )w1 Z (K+γ)w2 K  0 −a2 − (β Z+γ + 0 2 ¯ K+γ  1 1) ¯2 w w3 Z 3 ¯ 2 Z(c − ) 0 D D2 ¯ ¯ (β1 Z+γ (K+γ)w2 K 1 )w1 Z ¯ − w3 ) −a1 , λ2 = −a2 − (β + , λ3 = 2Z(c 2 ¯ K+γ D 1 Z+γ1 )

  

Neutral

−a1

31

Neutral =0

Table 8: Blow-up phenomena of system (2)-(6) Model c> (2)

w3 D

Trivial

Blow-up condition for Z w3 (non-delayed & delayed both) c< D 

|Z0 | ln |Z0 | ln



|Y0 | w3 −D c−δ1 |Y0 | w3 −D c−δ1

>



a2 +

w2 β

δ1

, |Y0 | >

w3 c−δ1

−D

, |Y0 | >

w3 c−δ1

−D

w

a2 + β 2

(3)

Trivial

(4)

non-delayed:Figure3(a), delayed: Figure3(b)   w2 a2 + w3 0| −D Trivial |Z0 | ln w|Y > δ1 β , |Y0 | > c−δ 3 −D 1

>

1

δ1

c−δ1

(5)

non-delayed:Figure4(a), delayed: Figure4(b)   w2 a + 2 w3 0| Trivial |Z0 | ln w|Y > δ1β1 , |Y0 | > c−δ −D 3 −D 1 c−δ1

(6)

non-delayed:Figure5(a), delayed: Figure5(b)   w2 a2 + w3 0| Trivial |Z0 | ln w|Y > δ1β1 , |Y0 | > c−δ −D 3 −D 1 c−δ1

non-delayed:Figure6(a), delayed: Figure6(b)

Table 9: Parametric values used for stability analysis of model system (2)-(6) Model (2)

(3)

(4)

(5)

(6)

Parameter set a1 = 1, K = 30, w = 0.75, b = 60, a2 = 0.05, w1 = 0.6, w2 = 0.2, w3 = 0.65, β = 0.21, γ = 50, c = 0.03, D = 1 a1 = 1, K = 30, w = 0.75, b = 60, a2 = 0.05, w1 = 0.6, w2 = 0.2, w3 = 0.65, β1 = 0.21, γ1 = 50, δ = 0.3, ρ = 0.03, c = 0.03, D = 1 a1 = 1, K = 30, w = 0.75, b = 60, a2 = 0.05, w1 = 0.6, w2 = 0.2, w3 = 0.5, β = 0.21, γ = 50, c = 0.01, D = 1, j = 16 a1 = 1, K = 30, w = 0.75, b = 60, a2 = 0.05, w1 = 0.6, w2 = 0.2, w3 = 0.65, β1 = 0.21, j = 18, γ1 = 50, δ = 0.3, ρ = 0.03, c = 0.03, D = 1 a1 = 1, K = 30, w = 0.75, β = 60, a2 = 0.05, w1 = 0.6, γ = 50, w2 = 0.2, w3 = 0.65, β1 = 0.21, γ1 = 50, δ = 0.3, ρ = 0.03, c = 0.03, D = 1

32

Critical value of τ τ0 = 3.4

Phase space diagram Figure 1(a)(i)

Bifurcation diagram Figure 1(a)(ii)

τ0 = 3.61

Figure 1(b)(i)

Figure 1(b)(ii)

τ0 = 11.617

Figure 1(c)(i)

Figure 1(c)(ii)

τ0 = 11.97

Figure 1(d)(i)

Figure 1(d)(ii)

τ0 = 3.45

Figure 1(e)(i)

Figure 1(e)(ii)

Table 10: Blow-up phenomena with parametric values of model system (2)-(6) Model (2)

(3)

(4)

(5)

(6)

Parameter set a1 = 1, K = 100, w = 1, b = 1.01, a2 = 0.01, w1 = 0.01, w2 = 10, w3 = 0.9, β = 10, γ = 0.2, c = 0.4, D = 0.5 a1 = 1, K = 100, w = 1, b = 1.01, a2 = 0.01, w1 = 0.01, w2 = 10, w3 = 0.9, β1 = 10, γ1 = 0.2, δ = 1, ρ = 0.2, c = 0.4, D = 0.5 a1 = 1, K = 100, w = 1, b = 1.01, a2 = 0.01, w1 = 0.01, w2 = 0.2, w3 = 0.9, β = 10, γ = 0.2, c = 0.4011, D = 0.5, j = .502 a1 = 1, K = 100, w = 1, b = 1.01, a2 = 0.01, w1 = 0.01, w2 = 10, w3 = 0.9, β1 = 10, j = .502, γ1 = 0.2, δ = 0.2, ρ = 0.2, c = 0.4011, D = 0.5 a1 = 1, K = 100, w = 1, β = 1.01, a2 = 0.01, w1 = 0.01, γ = 50, w2 = 10, w3 = 0.9, β1 = 10, γ1 = 0.2, δ = 0.2, ρ = 0.2, c = 0.40031, D = 1

Initial Condition (IC) (5, 6, 5)

Blow-Up without Delay No

Time series without Delay Figure 2(a)

Blow-Up with Delay No τ = 0.231

Time series with Delay Figure 2(b)

(5, 6, 5)

Yes

Figure 3(a)

Yes τ = 0.0231

Figure 3(b)

(5, 6, 5)

No

Figure 4(a)

No τ = 0.231

Figure 4(b)

(5, 6, 5)

Yes

Figure 5(a)

No τ = 1.9231

Figure 5(b)

(5, 6, 5)

Yes

Figure 6(a)

No τ = 0.0231

Figure 6(b)

33

(i)

(ii)

350

(i)

(ii)

450

300

400

140

200 350

120

250

τ=0 τ=2.5 τ=4

τ=0 τ=2.5 τ=4

150

60 40

200

300

Population

Z

Population

80

Z

100

100

150

250

200

50

20

E*

•

0 60

100

15

40

50

5 0

100 15

40

Z

10 20 Y

0

0

Y

X 0

1

2

τ0= 3.4

τ

4

20

5

5 0

X

0

0

(a) (ii)

120

35

15

60

•

Population

20

Z

25

τ

τ0= 3.61

4.3

(ii)

80

30 Population

Z

τ=0 τ=3 τ=12.4

40

80

35 E*

2

100

50 45

•

1

120

τ=0 τ=2.5 τ=12.2

30

0

(i)

100

50

40

Y X

(b)

(i)

45

Z

50

10

Y

X

150

E*

•

0 60

25 20 15

Y

60 Y

10

10 40

5

40

5

Z

0

Z

0 150

150 X

20

30

100

30

100

Y

50

10 0

X

20

20

20 Y

X 0

0

0

5

τ0=11.617

τ

15

50

10 0

X 0

0

(c)

0

5

τ

(d) (i)

(ii)

400

350 300 300

250

τ=0 τ=2.5 τ=4 250

Population

Z

200 150 100

200

150

50 0

100

•

60

15

40 Y

5 0

Y

50

10 20

Z

X 0

0

X 0

1

2 τ

τ0=3.45

4

(e)

Figure 1: (i) Phase space diagram of system (2)-(6) in (a,b,c,d,e) respectively, (ii) bifurcation diagram of system (2)-(6) in (a,b,c,d,e) respectively with the varying parameter τ . Other parameters are as in the table 9. 34

τ0=11.97

15

Model 1, No delay −I.C=(5,6,5) 160

120

120

100

100

80 60

80

60

40

40

20

20

0

0

0

500

1000

1500

2000

X(t) Y(t) Z(t)

140

Population density

Population density

140

−20

With delay-I.C=(5,6,5), lag=[0.231]

160

X(t) Y(t) Z(t)

-20

2500

0

50

100

150

200

250

300

350

400

time t

time t

(a)

(b)

Figure 2: (a) Blow-up does not occur in system (2) for same parameter set and IC as it does in Model (3) without delay, (b) blow-up does not occur in Model (2) for same parameter set and IC as it does in Model (3) with delay (τ = 0.231). Other parameters are as in the Table 10.

No delay −I.C=(5,6,5) 16000

3.5

X(t) Y(t) Z(t)

14000

X(t) Y(t) Z(t)

3

2.5

10000

Population density

Population density

12000

8000 6000

2

1.5

1

4000

0.5

2000 0

With delay-I.C=(5,6,5), lag=[0.0231]

×10 14

0

0.2

0.4

0.6

0.8

1

1.2

0

1.4

0

0.1

0.2

0.3

0.4

0.5

0.6

time t

time t

(a)

(b)

Figure 3: (a) Non delayed counter part of system (3) blows-up in finite time, (b) the delayed model also blows up with delay (τ = 0.231). Other parameters are as in the Table 10.

35

Model 3,No delay −I.C=(5,6,5) 160

120

120

100

100

80 60

80

60

40

40

20

20

0

0

0

50

100

150

200 time t

250

300

350

X(t) Y(t) Z(t)

140

Population density

Population density

140

−20

Model 3,With delay-I.C=(5,6,5), lag=[0.231]

160

X(t) Y(t) Z(t)

-20

400

0

500

1000

1500

2000

2500

3000

3500

4000

time t

(a)

(b)

Figure 4: (a) Blow-up does not occur in Model (4) for same parameter set and IC as it does in Model (5). This is without delay, (b) blow-up does not occur in Model (4) for same parameter set and IC as it does in Model (5) with delay (τ = 0.231). Other parameters are as in the Table 10.

Model 4 ,No delay −I.C=(5,6 ,5) 18000

14000

14000

12000

12000

10000 8000

10000

8000

6000

6000

4000

4000

2000

2000

0

0.2

0.4

0.6

0.8

1

1.2

X(t) Y(t) Z(t)

16000

Population density

Population density

16000

0

Model 4 With delay-I.C=(5,6,5), lag=[1.9231]

18000

X(t) Y(t) Z(t)

0

1.4

0

0.5

1

1.5

2

2.5

3

time t

time t

(a)

(b)

Figure 5: (a) Non-delayed system for Model (5) blows-up, (b) the short time dynamics with delay (τ = 1.9231) for Model (5), blow-up is avoided, (c) the long time dynamics with delay (τ = 1.9231) for Model (5), blow-up is avoided and the delay is providing a damping effect. Other parameters are as in the Table 10.

36

Model 5,No delay −I.C=(5,6,5)

4

2.5

x 10

X(t) Y(t) Z(t)

5000

2

4000 Population density

Population density

With delay-I.C=(5,6,5), lag=[0.0231]

6000

X(t) Y(t) Z(t)

1.5

1

3000

2000

0.5

0

1000

0

0.2

0.4

0.6

0.8

1

1.2

0

1.4

0

0.5

1

1.5

2

2.5

3

time t

time t

(a)

(b)

Figure 6: (a) Non-delayed system for system (6) blows-up, (b) the short time dynamics with delay (τ = 0.231) for system (6), blow-up is avoided and the delay is providing a damping effect and (c) the long time dynamics with delay for Model (6), notice how the delay is providing a damping effect. Others parameters are as in the Table 10.

37

30

30

25

25

20

20

15

15

AS

10

5

BL

Z0

Z0

BL

1

AS

10

2

3

4

5

5

1

2

(a)

3

4

5

4

5

(b)

0.8

30

25 0.6

20

BL

Z0

Z0

BL 0.4

15

0.2

0

AS

1

AS

10

2

3

Y

4

5

5

1

2

3

Y

0

(c)

0

(d)

20

16

BL

Z0

12

8

AS

4

0

1

2

3

4

5

Y0

(e)

Figure 7: Stability region with varying initial values Y0 and Z0 of model systems (2)-(6). Here ’AS’ and ’BL’ represent Asymptotically stable and Blow up respectively. Others parameters are as in the Table 10. 38