Modeling the Allee effects induced by cost of predation fear and its carry-over effects
MModeling the Allee effects induced by cost of predation fear and itscarry-over effects
Sourav Kumar Sasmal a, ∗ , Yasuhiro Takeuchi b a Department of Mathematics, Birla Institute of Technology and Science, Pilani, Rajasthan 333-031, India b Department of Physics and Mathematics, Aoyama Gakuin University, Kanagawa 252-5258, Japan
Abstract
Predation driven Allee effects play an important role in the dynamics of small population, however,such predation-driven Allee effects cannot occur for the model with type I functional response. Itgenerally occurs when a generalist predator targets some specific prey. However, apart from the lethaleffects of predation, there are some non-lethal effects in the presence of predator. Due to the fear ofpredation, positive density dependence growth may be observed at low population density, because ofreduced foraging activities. Moreover, this non-lethal effect can be carried over generations. In thepresent manuscript, we investigate the role of predation fear and its carry-over effects in prey-predatormodel. First, we study the single species model in global perspective. We have shown that dependingon the birth rate, our single species model describes three types of growth dynamics, namely, strongAllee dynamics, weak Allee dynamics and logistic dynamics. Then we consider the explicit dynamics ofpredator, with type I functional response. Basic dynamical properties, as well as global stability of eachequilibria have been discussed. From our analysis, we can observe that both the fear and its carry-overeffects have significant role in the stability of the coexistence equilibrium, even if for the model withtype I functional response. The phenomenon ’paradox of enrichment’ can be observed in our model,which cannot be observed in the classical prey-predator model with type I functional response. However,we can see that such phenomenon can be ruled out by choosing suitable non-lethal effect parameters.Therefore, our study shows how non-lethal effects change the dynamics of a competition model, and hasimportant biological insights, specially for the understanding of the dynamics of small populations. Keywords:
Allee effect, fear effect, carry-over effect, predator-prey interactions, global stability.
1. Introduction
Allee effect, is a positive density-dependence phenomena, which is defined as the positive relation-ship between population density and per capita growth rate ( pgr ) at low population density (Allee, 1931;Courchamp et al., 2008). In contrast to the logistic dynamics, Allee effects play an important role tothe extinction of small populations. There are a number of mechanisms for which Allee effects havebeen observed, such as, mate limitation, cooperative defense, cooperative feeding, environmental condi-tioning, inbreeding depression, demographic stochasticity, etc. (Courchamp et al., 2008; Stephens et al.,1999; Kramer et al., 2009; Dennis, 1989; Lewis and Kareiva, 1993). Apart from the above mechanisms,predator-driven Allee effects also can be observed in nature. Though, predator-driven Allee effects are ∗ Corresponding author.
Email addresses: [email protected], [email protected] (Sourav Kumar Sasmal), [email protected] (Yasuhiro Takeuchi)
Preprint submitted to arXiv January 27, 2021 a r X i v : . [ q - b i o . P E ] J a n imited as predator population declines, prey population also declines - leading to the negative densitydependence growth. It can be observed for some specific type of functional response, such as Hollingtype II functional response (Courchamp et al., 2008; Gascoigne and Lipccius, 2004; Kramer and Drake,2010). However, here we will investigate the occurrence of Allee effects for Holling type I functionalresponse, in the presence of non-lethal effects of predation. Due to considerable impact of Allee effects,it seeks significant attention among theoretical ecologists in several aspects, like population ecology(Courchamp et al., 2008; Dennis, 1989, 2002; Sasmal and Ghosh, 2017; Sasmal et al., 2017), biologicalinvasion (Drake, 2004; Lewis and Kareiva, 1993; Taylor and Hastings, 2005), eco-epidemiology (Deredecand Courchamp, 2006; Kang et al., 2014; Sasmal and Chattopadhyay, 2013) etc.Apart from the direct killing, which is widely observed in nature, many preys modify their traitsin response to the predation risk. These modified traits could be related to behavior, morphology, lifehistory of prey. To avoid predation, prey shows various anti-predator behaviors, e.g., habitat changes,reduced foraging activities, vigilance, some physiological changes, etc. (Cresswell, 2011). Such effectsare known as trait-mediated indirect effect, as such effect arises from predator’s influence on prey traits,rather than arises from prey densities. Such non-lethal predator effects could be immediate and caninfluence entire prey population over entire lifetime. It can be argued that non-lethal effects are impor-tant and are needed to consider in population ecology only when it is large compared to direct densitydependent effects. However, it may not be the case even when predation rate is high. Many authorssuggest that such indirect effects could be equivalent or more influential compared to density effects(direct predation) (Creel and Christianson, 2008; Cresswell, 2011; Preisser and Bolnick, 2008). Thisargument was supported by the experimental data from prey-predator interaction of larval dragonfly - Anax sp. (predator) and bullfrog tadpoles -
Rana catesbeiana (prey) (Peacor and Werner, 2001). Somerecent studies showed that among such anti-predator responses, fear of predation can play an importantrole as direct predation effect in prey-predator models (Zanette et al., 2011; Wang et al., 2016; Sasmal,2018). Due to the predation fear, scared prey forages less, as well as it leads to some stress-relatedphysiological changes, which impact reproduction success (Creel and Christianson, 2008; Schmitz et al.,1997; Elliott et al., 2016). For example, birds flee from their nests in response to their predators sound asan anti-predator response (Creel and Christianson, 2008; Cresswell, 2011). Though such anti-predatorresponse may be instantly beneficial as it increases the adult survival, however, as a long-term cost itreduces the reproduction rate (Cresswell, 2011). Zanettee et al. (Zanette et al., 2011) experimentallyshowed the 40% reduction in offspring production of prey (song sparrows -
Melospiza melodia ) due topredation fear by providing predatory sound only and without direct killing. They showed that thisreduction is due to anti-predator behavior which affects the reproduction of song sparrows. Therefore,for free living wildlife population, incorporation of non-lethal trait-mediated indirect effect by predationfear is important. Moreover, as density declines, prey individuals are more vigilant and less foraging.Therefore, such non-lethal effect could be a cause of Allee effects and increases the extinction risk ofsmall populations (Clutton-Brock et al., 1999; Mooring et al., 2004).The ’carry-over effect’ originally started from recurrent measures of clinical experiments. However,recently it has been used in ecological and evolutionary aspects and can be used for broad range ofsituations. O’Connor et al. (O’Connor et al., 2014) proposed the following working definition for carry-over effects: “In an ecological context, carry-over effects occur in any situation in which an individual’sprevious history and experience explains their current performance in a given situation”. In view of theabove definition, carry-over effects are not restricted to the seasonal requirement, discrete-time scale,migration, etc. (Marshall and Morgan, 2011; O’Connor et al., 2011), which was previously considered.2t should be considered as a more general phenomenon, which allows us to identify in broad range ofsituations, like, within and across life-history stages, seasons, years etc. Under this definition, life-historytrade-offs and costs of reproduction can be viewed as special types of carry-over effects. Moreover, somelab experiments showed that non-lethal carry-over effects have impact in long-term population dynamics(Betini et al., 2013a,b). Carry-over effects can occur in both across multiple seasons or within a singleseason (e.g., transitions between physiological states within a season). Experimental evidences of carry-over effects within a single season and over short time periods are observed in insects (De Block andStoks, 2005), amphibians (Touchon et al., 2013), marine fish (Johnson, 2008), and marine invertebrates(Marshall and Morgan, 2011), etc. Due to the above reasons, study on ecological carry-over effects hasbeen increasing in mathematical modeling studies (Norris, 2005; Runge and Marra, 2005; Norris andTaylor, 2006), as well as in empirical research (Norris et al., 2004; Inger et al., 2010; Legagneux et al.,2012; Sedinger et al., 2011). Therefore, integrating the research on carry-over effects, with potentialconnection between life-history trade-offs and cost of reproduction will improve our understanding ofthe factors affecting population dynamics in nature.Betini et al. (Betini et al., 2013a) introduced an experimental model system of
Drosophila to studythe sequential density dependence and carry-over effects and used a simple
Ricker map with season-specific parameters. On the other hand, Elliott et al. Elliott et al. (2017) investigated the role offear in relation to fitness and population density by considering a prey-predator system of
Drosophilamelanogaster (prey) and mantid (predator), both in breeding and non-breeding seasons. They used theexperimental results to parameterize a bi-seasonal Ricker map and provided the evidence that indirecteffect of predator can be a cause of Allee effect, which is very important to understand the dynamicsof small population. Motivated from the above discussion, we develop and analyze continuous-timepopulation models to investigate the cost of predation fear and its carry-over effects in prey-predatorinteraction with Holling type I functional response. The objectives of our study are to answer thefollowing questions: ( i ) How these non-lethal effects (cost of predation fear and its carry-over effects)change the growth dynamics of prey population? ( ii ) How the phenomenon ’paradox of enrichment’ can be observed for our model and under which condition it can be ruled out? ( iii ) Role of non-lethaleffects on the stability of the coexistence equilibrium. ( iv ) What is the global dynamical behavior ofour proposed model? The remainder of the paper is organized as follows: In Section 2, we formulateand analyze a single species population model with fear and its carry-over effects. We find the globaldynamical behavior of our proposed single species model, as well as we show that how growth dynamicschanges depending on the parameter values. In Section 3, we consider the explicit dynamics of predatorpopulation, with Holling type I functional response. Basic dynamical properties, as well as globalstability of equilibria are provided in this section. Moreover, existence and stability of Hopf-bifurcationis also provided in this section. Some numerical simulation results are shown in the Section 4. In Section5, we discuss our results and findings and provide some potential future directions. Some detailed proofsof our analytical findings are given in Section 6.
2. Single species model with fear and carry-over effects
First, we consider prey growth follows the logistic dynamics, which can be split into three parts,birth, natural death, and death due to intra-prey competition. Thus in the absence of predator, a singlespecies population model is given by the following ODE: dxdt = rx − d x − d x , (2.1)3here x is the density of prey, whose maximum birth rate is r in the absence of predation. d is thenatural death rate and d is the density dependent death rate of prey. Next, we consider a single species(prey) population model with a generalist predator, which is at constant density y . Here, we neglect thedirect predation. We develop and analyze a single species population model with fear and carry-overeffects which is given by dxdt = rx (cid:124)(cid:123)(cid:122)(cid:125) birth cx cx + f y (cid:124) (cid:123)(cid:122) (cid:125) fear and carry-over effect − d x − d x = x (cid:20) r (1 + cx )1 + cx + f y − d − d x (cid:21)(cid:124) (cid:123)(cid:122) (cid:125) per capita growth rate ≡ x Θ( x ) . (2.2)Here c is the carry-over effect parameter due to fear, which is quantified by the parameter f . In theModel (2.2), if f = 0 (i.e., if there is no growth-rate reduction due to predation fear), then the modelwill be simply logistic dynamics (2.1). If c = 0, then model is reduced to the single species model withonly fear effect which has been studied in (Wang et al., 2011; Sasmal and Takeuchi, 2020). Lemma 2.1.
The solution of system (2.2) is uniformly ultimately bounded in R + with lim t →∞ x ( t ) ≤ r − d d , when r ≥ d . Here, System (2.2) can exhibit Allee effects depending on the parameters r , c , f , d , d , and whenΘ (cid:48) (0) >
0. The Allee effect will be weak if population of (2.2) persists with Θ(0) > <
0. Moreover, when Θ (cid:48) (0) < f and c is zero (or both) then the system (2.2) shows logistic dynamics.System (2.2) always has the extinction equilibrium x = 0. Other equilibria of the System (2.2) areroots of the quadratic equationΦ( x ) ≡ cd x + [ c ( d − r ) + d (1 + f y )] x + [ d − r + d f y ] = 0 . (2.3)We rename the coefficients of the above equation asΩ = c ( d − r ) + d (1 + f y ) and Ω = d − r + d f y. Denote the roots of the above equation by x , = − Ω ± √ ∆2 cd , where ∆ = Ω − cd Ω = [ c ( d − r ) + d (1 + f y )] − cd [ d − r + d f y ] . (2.4)Here, x > x (when both the roots are real, i.e., when ∆ > x ) < x >
0, when r ≤ d . Hence, x ( t ) → t → ∞ , when r ≤ d . No other non-negative equilibrium ex-ists in this case, except x = 0, which is globally asymptotically stable. Hereafter, we assume that r > d .Depending on the number of roots of the quadratic equation (2.3), we have the following three cases:4. (2.3) has no positive solution in the following two cases:(i) Ω ≥ ≥ , (ii) Ω < <
02. (2.3) has unique positive solution in the following three cases:(i) Ω < , (ii) Ω < = 0 , (iii) Ω < .
3. (2.3) has two positive solutions if Ω < , Ω > > . Following theorem summarizes the existence and stability of all equilibria for the Model (2.2).
Theorem 2.1. [Existence and stability of equilibria for Model (2.2) .]1. The System (2.2) has only the trivial extinction equilibrium x = 0 , if any of the following twoconditions holds:(i) r ≤ min (cid:110) d (1 + f y ) , d + d (1+ fy ) c (cid:111) (see Figure 1(a)).(ii) r > d + d (1+ fy ) c and ∆ < (see Figure 1(b)).In this case the equilibrium x = 0 is globally asymptotically stable.2. The System (2.2) has two non-negative equilibria, namely, the trivial extinction equilibrium x = 0 ,and a non-trivial equilibrium, if any of the following three conditions holds:(i) If r > d (1 + f y ) , then (2.2) has unique positive equilibrium point x . Here, x is alwaysunstable and x is globally asymptotically stable (see Figure 1(c)).(ii) If r = d (1 + f y ) > d + d (1+ fy ) c , then (2.2) has unique positive equilibrium x = c ( r − d ) − d (1+ fy ) cd . Here also x is unstable and x is globally asymptotically stable (see Figure1(d)).(iii) If r > d + d (1+ fy ) c and ∆ = 0 , then (2.2) has one positive equilibrium x = c ( r − d ) − d (1+ fy )2 cd of order two. Here, x is locally asymptotically stable and x is a saddle (see Figure 1(e)).3. The System (2.2) has three non-negative equilibria, namely, the trivial extinction equilibrium x ,and two non-trivial equilibria x and x , if d + d (1+ fy ) c < r < d (1 + f y ) and ∆ > . Here, boththe equilibria x and x are locally asymptotically stable, and x is always unstable (see Figure1(f )).Where ∆ is defined in (2.4) . Note:
Under conditions of Theorem (2.1) part 1( i ) and 1( ii ), pgr (per capita growth rate) of theModel (2.2) is always negative, (Θ( x ) < x > xdx/dtx (a) No positive equilibrium, x is globally stable for r = 1 .
25. Here, d (1 + fy ) = 2 and d + d (1+ fy ) c =1 . xdx/dtx (b) No positive equilibrium, x is globally stable for r = 1 .
85. Here, d + d (1+ fy ) c = 1 . − . < dx/dtx x x (c) Unique positive equilibrium x is globally stablefor r = 2 .
25. Here, d (1 + fy ) = 2. dx/dt xx x (d) Unique positive equilibrium x is globally stablefor r = d (1 + fy ) = 2. Here, d + d (1+ fy ) c = 1 . -3 dx/dt xx x (e) Unique positive equilibrium x is unstable and x is locally stable for r = 1 . d + d (1+ fy ) c =1 .
5, and ∆ ≈ dx/dt xx x x (f) Two positive equilibria x and x exist for r = 1 . x and x are locally stable. Here, d (1 + fy ) = 2, d + d (1+ fy ) c = 1 .
5, and ∆ = 0 .
06 ( > • ) are globally stable equilibria, reddots ( • ) are unstable equilibria, yellow dots ( • ) are locally stable equilibria, and blue arrows ( → ) are the flow direction.The fixed parameters are f = 1, c = 1, y = 1, d = 1, d = 0 . igure 2: Bifurcation diagram of our Model (2.2), with respect to the parameter r . Saddle-node bifurcation occurs at r = 1 .
866 (black dot-dashed line) and transcritical bifurcation occurs at r = 2 (black dashed line). All the other parametersare fixed at f = 1, y = 1, c = 1, d = 1, and d = 0 . Type of growth dynamics Parameter constraints
Weak Allee dynamics r > max (cid:110) d (1 + f y ) , d (1+ fy ) cfy (cid:111) Strong Allee dynamics max (cid:110) d + d (1+ fy ) c , d (1+ fy ) cfy (cid:111) < r < d (1 + f y ) and ∆ > r < d (1+ fy ) cfy or c = 0 or f = 0 Table 1: Type of growth dynamics and parameter constraints for the Model (2.2).
Lemma 2.2. [Saddle-node and transcritical bifurcations.] When the Model (2.2) parameters are suchthat ∆ = 0 , then a saddle-node bifurcation occurs at the positive equilibrium point x . Moreover, when Ω = 0 , i.e., r = d (1 + f y ) , then Model (2.2) experiences a transcritical bifurcation at the positiveequilibrium x . The saddle-node bifurcation and transcritical bifurcation have shown in the Figure 2.For the Model (2.2), the pgr is given by Θ( x ) and Θ (cid:48) ( x ) = rcfy (1+ cx + fy ) − d . Therefore, accordingto the previous discussion, Model (2.2) may show Allee dynamics if Θ (cid:48) (0) >
0, i.e., if r > d (1+ fy ) cfy .Therefore, our model (2.2) shows three types of growth dynamics, which is summarized in the Table (1). • If r > d (1 + f y ), then the population of species is persistent in R + . Moreover, Θ(0) ≥ r ≥ d (1 + f y ). Therefore, the Model (2.2) shows weak Allee dynamics (Figure 3(a)) if r > max (cid:110) d (1 + f y ) , d (1 + f y ) cf y (cid:111) . • If d + d (1+ fy ) c < r < d (1 + f y ) and ∆ >
0, then the population of species is persistent in R + \ [0 , x ), and goes to extinction otherwise. Since Θ(0) < r < d (1 + f y ). TheModel (2.2) shows strong Allee dynamics (Figure 3(b)) if the following conditions holdsmax (cid:110) d + d (1 + f y ) c , d (1 + f y ) cf y (cid:111) < r < d (1 + f y ) and ∆ > . x is known as Allee threshold. Population will persist if initial population density is above x , otherwise it goes extinction. • Moreover, if r < d (1+ fy ) cfy , then the System (2.2) shows logistic growth dynamics (Figure 3(c)).Also, if f = 0 or c = 0 (or both are zero), then also equation (2.2) shows logistic growth dynamics(Figure 3(d)). Remark 2.1.
In this section, we do not consider the predator dynamics explicitly. We only consider theconstant predator population, without predation term. However, the qualitative properties of the Model (2.2) , will be the same for constant predator population model with Holling type I functional response.In fact, we rewrite the model as, dxdt = rx cx cx + fy − d (cid:124)(cid:123)(cid:122)(cid:125) constant death x − d x − axy (cid:124)(cid:123)(cid:122)(cid:125) type- I response = rx cx cx + fy − ( d + ay ) (cid:124) (cid:123)(cid:122) (cid:125) constant death x − d x . In the next section, we will discuss the dynamics in the presence of predator population explicitly,with Holling type- I functional response.
3. Predator-prey model with type - I response function
In this section, we study a prey-predator model with fear and carry-over effects with linear functionalresponse (Holling type- I ). Thus, our two-species prey-predator model becomes, dxdt = rx (1+ cx )1+ cx + fy − d x − d x − axy dydt = aαxy − my. (3.1)Here, a is the rate of predation, α is the conversion efficiency from prey biomass to predator biomassand m is the natural death rate of predator population. In the following theorem we summarize thebasic dynamical properties of the Model (3.1). Lemma 3.1. [Positivity and boundedness of solutions for Model (3.1) .] For the System (3.1) , the set R is positively invariant. Moreover, the System (3.1) is dissipative, i.e., every solution of (3.1) isultimately bounded in R , with the following properties lim t →∞ sup x ( t ) ≤ r − d d lim t →∞ sup (cid:2) x ( t ) + α y ( t ) (cid:3) ≤ (cid:26) r − d d if m>r − d r − d m )24 d m if m ≤ r − d . The Model (3.1) always has the trivial extinction equilibrium E = (0 , E = (cid:16) r − d d , (cid:17) , under the condition r > d . Moreover,there exists a unique interior equilibrium E ∗ = ( x ∗ , y ∗ ) = (cid:0) maα , y ∗ (cid:1) where y ∗ is the root of the quadraticequation Γ y + Γ y + Γ = 0 , Density (x) -0.100.10.2
PGR ((x)) (a) Model (2.2) shows weak Allee dynamics when r (=2 . > max { d (1 + fy ) , d (1+ fy ) cfy } (= max { , } ) forthe parameter values r = 2 . f = 1 and c = 1. Density (x) -0.2-0.100.1
PGR ((x)) (b) Model (2.2) shows strong Allee dynamics when d (1+ fy ) cfy (= 1) < r (= 1 . < d (1 + fy )(= 2) forthe parameter values r = 1 . f = 1 and c = 1. Density (x) -0.200.20.40.6
PGR ((x)) (c) Model (2.2) shows logistic dynamics when d (1 + fy )(= 1 . < r (= 2) < d (1+ fy ) cfy (= 4 .
69) for theparameter values r = 2, f = 0 . c = 0 . Density (x) -0.500.51
PGR ((x)) (d) Model (2.2) shows logistic dynamics when either f = 0 or c = 0 (or both are zero) for the parametervalues r = 2, f = 0 . c = 0 .
3. Solid line is for f = 0 (or when both are zero) and dot-dashed line isfor c = 0.Figure 3: Other parameters are fixed as y = 1, d = 1 and d = 0 .
25. Dashed line represent the x − axis (i.e., when PGRis zero). = a α f Γ = aα [ a ( aα + cm ) + f ( d aα + d m )]Γ = ( aα + cm )( d aα + d m − raα ) . The above equation has unique real positive root y ∗ = − [ a ( aα + cm )+ f ( d aα + d m )]+ √ [ a ( aα + cm ) − f ( d aα + d m )] +4 a αrf ( aα + cm )2 a αf , (3.2)iff raα > d aα + d m i.e., iff r > d + d maα . (3.3) Remark 3.1.
When the prey and predator population coexists at E ∗ , the density of prey population doesnot depend on any of the parameters, cost of predation fear ( f ) and carry over effect ( c ) . However, thedensity of predator population depends on both the parameters. We have, ∂y ∗ ∂f = − ( aα + cm ) (cid:104) a ( aα + cm ) − f ( d aα + d m )+2 aαrf − √ [ a ( aα + cm ) − f ( d aα + d m )] +4 a αrf ( aα + cm ) (cid:105) aαf √ [ a ( aα + cm ) − f ( d aα + d m )] +4 a αrf ( aα + cm ) . (3.4) Now, it is easy to prove that the numerator of (3.4) is positive, iff r > d + d maα , which is the existencecondition of y ∗ . Thus, ∂y ∗ ∂f < , i.e., at the coexistence state, as we increase the cost of fear, the densityof predator population decreases.Similarly, we have ∂y ∗ ∂c = − m (cid:104) √ [ a ( aα + cm ) − f ( d aα + d m )] +4 a αrf ( aα + cm ) − a ( aα + cm )+ f ( d aα + d m ) − aαrf (cid:105) a √ [ a ( aα + cm ) − f ( d aα + d m )] +4 a αrf ( aα + cm ) . (3.5) Again, it is easy to prove that the numerator of (3.5) is negative, iff r > d + d maα , which is the existencecondition of y ∗ . Thus, ∂y ∗ ∂c > , i.e., at the coexistence state, as we increase the carry over effect, thedensity of predator population increases. The following theorem describes the local stability of all three equilibria.
Theorem 3.1. [Local stability of equilibria for Model (3.1) .]1. The extinction equilibrium E is locally asymptotically stable if r < d , and a saddle otherwise.2. The prey-only equilibrium E is locally asymptotically stable if the condition (3.3) is reversed anda saddle when the condition (3.3) is satisfied.3. The coexistence equilibrium E ∗ is locally asymptotically stable if r < d ( aα + cm + aαfy ∗ ) a α cfy ∗ , where y ∗ is defined in (3.2) .Moreover, if d + d maα < r < d ( aα + cm ) aαc , then the coexistence equilibrium E ∗ is always locallyasymptotically stable. The following theorems give the additional stability properties for the coexistence equilibrium.10 heorem 3.2. If r < r c , then the coexistence equilibrium E ∗ is always locally asymptotically stable.Furthermore, if r > r c , then the sufficient condition for the local stability of the coexistence equilibrium E ∗ is f < f c , where r c and f c are given by r c = c ( aα + cm )( d aα + d m )+ √ c ( aα + cm )( d aα + d m )[ c ( aα + cm )( d aα + d m )+4 a α d ]2 aαc ( aα + cm ) ,f c = a αd [ aαd + rc ( aα + cm )] rc [ aαr c ( aα + cm ) − rc ( aα + cm )( d aα + d m ) − aαd ( d aα + d m )] . Remark 3.2. If r < r c , then the condition (6.2) is always satisfied, i.e., when r is small (smaller than r c ), the fear parameter f has no role in the coexistence equilibrium stability. In other words, if d + d maα < r < r c , then the equilibrium E ∗ is always locally asymptotically stable. The stability of E ∗ will not change if thebirth rate of prey is not large enough to support oscillations, which is similar to the result obtained byWang et al. (2016) with type- II response function.If r > r c , then the condition (6.2) is satisfied if f < f c , which is the sufficient condition for thestability of E ∗ , when r is large (larger than r c ). In other words, if r > max (cid:110) d + d maα , r c (cid:111) , then the sufficient condition for the local stability of E ∗ is f < f c .Therefore, we can see that the fear parameter has an important role on the stability of the coexistenceequilibrium, even for the type- I response function. However, the actual stability property with respect tothe fear parameter f is more complex than the above scenarios, which is discussed later through numericalsimulations. Theorem 3.3. If r < d maα , then the coexistence equilibrium E ∗ is locally asymptotically stable. Fur-thermore, if r > d maα , then the sufficient condition for the local stability of E ∗ is c < c r , which does notdepend on the fear parameter f , where (cid:16) c r = d aαraα − d m (cid:17) . Remark 3.3. If r < d maα , then the condition (6.5) is always satisfied, i.e., when r is small, theparameter c has no role in the stability of the coexistence equilibrium. In other words, if d + d maα < r < d maα , then the equilibrium E ∗ is always locally asymptotically stable. The stability of E ∗ will not be changed ifthe birth rate of prey is not large enough to support oscillations.If r > d maα , then the condition (6.5) is satisfied if c < c r , which is the sufficient condition for thestability of E ∗ . In other words, if r > max (cid:110) d + d maα , d maα (cid:111) , then the sufficient condition for the local stability of E ∗ is c < c r (6.6) . r is small,the stability of the coexistence equilibrium is not affected by the cost of fear or carry over effect. Forthe classical prey-predator model with Holling type I functional response, without cost of fear and carryover effect, Hopf-bifurcation never occur (our Model (3.1) will be reduced to the classical prey-predatormodel with Holling type I functional response, if we neglect the cost of fear). The result is same forthe prey-predator model with Holling type I functional response with only the cost of fear. However,for our model with both cost of fear and carry over effect, Hopf-bifurcation occurs as we increase theparameter r , and the phenomenon ’paradox of enrichment’ appears (McAllister et al., 1972; Riebesell,1974; Rosenzweig, 1971; Gilpin and Rosenzweig, 1972). When birth rate of prey is large enough, preyand predator can still go to the coexistence steady state according to preys cost of fear or carry overeffect. If either of the cost of fear and carry over effect is small enough, then it can suppress oscillations.Therefore, by incorporating the cost of fear and carryover effect, the phenomenon ’paradox of enrich-ment’ can occur, however, we can rule out such phenomenon by choosing suitable f or c . It is to benoted that the carrying capacity (e.g. for logistic model) is considered as the parameter to be evaluatedin terms of enrichment (Morozov et al., 2007). In this study, the parameter r , which is the maximumbirth rate of prey, is considered as the parameter to be evaluated in terms of enrichment. Actually, ifwe simplify the Model (2.1) then we will get the expression for the carrying capacity as r − d d . Fromthis expression, we can see that the carrying capacity increases or decreases with the parameter r , whenother parameters are fixed. Therefore, it is reasonable to assume the parameter r as the parameter tobe evaluated in terms of enrichment.All the local stability conditions in the Theorem (3.1) are actually global conditions. In the nexttheorems we will discuss about the global stability of all equilibria for Model (3.1). Theorem 3.4. [Global stability of boundary equilibria for Model (3.1) .] The equilibrium E is globallyasymptotically stable if r ∈ (0 , d ) and E is globally asymptotically stable if r ∈ (cid:0) d , d + d maα (cid:1) . Theorem 3.5. [Global stability of interior equilibrium for Model (3.1) .] The positive equilibrium E ∗ isglobally asymptotically stable if r ∈ (cid:0) d + d maα , d c (cid:1) . Remark 3.4.
From the Theorem (3.1) , we can see that, as the parameter r increases the Model (3.1) experiences two bifurcation of equilibrium and a Hopf-bifurcation at positive equilibrium (discussed later).When < r < d , E is globally asymptotically stable, when r passes d , E loses its stability to E ,which becomes globally asymptotically stable in d < r < d + d maα . Again, when r passes d + d maα then E loses its stability to E ∗ , which is locally asymptotically stable in d + d maα < r < d ( aα + cm + aαfy ∗ ) a α cfy ∗ , where y ∗ is defined in (3.2) . Moreover, when r passes through d ( aα + cm + aαfy ∗ ) a α cfy ∗ , then E ∗ loses its stability andlimit cycle oscillation occurs around E ∗ through Hopf-bifurcation.3.1. Existence of limit cycles and Hopf-bifurcation In this subsection, we investigate the possibility of Hopf-bifurcation at the coexistence equilibrium E ∗ by considering the fear effect parameter f , as the bifurcation parameter. Similarly, we can obtainthe Hopf-bifurcation with respect to other parameters r and c .At the Hopf-bifurcation point, the real parts of the eigenvalues of the characteristic equation (6.1)equal to zero. We set, at f = f H the Hopf-bifurcation occurs, which givesΨ ( f H ) = 0 and Ψ ( f H )Ψ ( f H ) > . ( f H ) = 0 ⇒ f H = aαrc − d ( aα + cm ) ± (cid:112) aαrc [ aαrc − d ( aα + cm )]2 d aαy ∗ . Moreover, if we simplify the above condition by using maple software, we can obtain that the valuesof f H are the roots of the quadratic equation C f + C f + C = 0 , where C = d [ c ( d aα + d m )( arα − ( d aα + d m )) − arαd ( aα + cm )] C = ac [ arcα ( arα − ( d aα + d m )) + d ( aα + cm )(2( d aα + d m ) − arα )] C = − a cd ( aα + cm ) . Therefore, f H = − C ± √ C − C C C . Further, at the Hopf-bifurcation point, we have d ( λ r ) df (cid:12)(cid:12)(cid:12) f = f H = 2 λ i d (Ψ ) df − Ψ d (Ψ Ψ ) df Ψ + 4 λ i (cid:12)(cid:12)(cid:12) f = f H (cid:54) = 0 , which is true if (cid:104) λ i d (Ψ ) df − Ψ d (Ψ Ψ ) df (cid:105) (cid:12)(cid:12)(cid:12) f = f H (cid:54) = 0, where λ r and λ i are real and imaginary parts ofthe eigenvalues of the characteristic equation (6.1).The following theorem gives conditions for the existence of Hopf-bifurcation at E ∗ for Model (3.1). Theorem 3.6. [Condition for the existence of Hopf-bifurcation at interior equilibrium for Model (3.1) .]If Ψ ( f H )Ψ ( f H ) > , and (cid:104) λ i d (Ψ ) df − Ψ d (Ψ Ψ ) df (cid:105) (cid:12)(cid:12)(cid:12) f = f H (cid:54) = 0 , hold, then the interior equilibrium E ∗ of Model (3.1) is locally asymptotically stable when f < f H , and undergoes Hopf-bifurcation at E ∗ when f = f H . The following theorem gives the direction and stability of Hopf-bifurcation around the coexistencesteady state E ∗ . Theorem 3.7. [Direction and stability of Hopf-bifurcation at interior equilibrium for Model (3.1) .] Letdefine L as L := 3 g v g uuu [(1 − g uv z ) + 2 z ( g uvv z w − g uuv z g v )] + 3 wg vvv z [2 w ( h u g uvv z − wg uuv z ) + g uv h u ]+ g uvv w [ w (1 − g uv z ) + 2 z ( g v g uu − h u g vv ) + 2 z w ( wg uvv z − g v z g uuv )]+ 2 g uuv z (cid:2) g v ( wg uuv z − g uu ) − w ( g vv + g v z g uvv ) (cid:3) + g uv ( h u g vv + wg v g uuv z − g v g uu ) . Then the Hopf-bifurcation is supercritical if
L < and it is subcritical if L > . r Prey (a) Bifurcation diagram of prey w.r.t. r . r Predator (b) Bifurcation diagram of predator w.r.t. r .Figure 4: Bifurcation diagrams for Model (3.1) with respect to r . Here, r varies from 0 to 2 .
5, with the initial condition[ x (0) , y (0)] = [0 . , . f = 1, c = 0 . d = 0 . d = 0 . a = 0 . α = 0 .
5, and m = 0 .
4. Numerical simulations
In the Figure (4), we fix the parameters f , c , d , d , a , α and m (the specific values are given in thefigure) and vary the parameter r . From the figure, we can see that as we increase the parameter r , thecoexistence equilibrium becomes unstable through Hopf-bifurcation and prey-predator oscillation occurs.This phenomenon is known as ’paradox of enrichment ,’ which can not be observed in the classical prey-predator model with Holling type I functional response (even in the presence of only fear effect (Wanget al., 2016; Sasmal, 2018)). Here, we consider the parameter r as the enrichment parameter because,if we simplify the Model (2.1), we will get the expression for carrying capacity as r − d d , therefore, it isreasonable to assume r as the enrichment parameter, when other parameters are fixed.In Theorems (3.2) and (3.3), we found the critical values of r and have shown that if r is less thanthe critical values then oscillation cannot be observed for our model (3.1). The coexistence equilibriumis always stable irrespective of the non-lethal parameters values. For oscillation, r must be greaterthan some threshold value, and in that case also, oscillation can be suppressed by the non-lethal effectsparameters, which is shown in Figures (5) and (6). Theoretically, we prove that for large values of r , if f or c is sufficiently small, then we can suppress the oscillation, however, the situation is morecomplex, which is shown in these two figures. In both the Figures (5) and (6), we choose r sufficientlylarge (greater than the critical values) and observed that predator-prey oscillation can occur only forthe intermediate values of f and c . Therefore, oscillation can be suppressed by both sufficiently low andhigh values of the cost of fear and carry over effect parameters, and we can rule out the phenomenon ’paradox of enrichment.’ We draw phase diagrams in Figure (7), to show how we can rule out theoscillating behavior by choosing high anti-predator response. Moreover, in Figure (8), we have shownthe Hopf-bifurcation curve for our Model (3.1), in the r − f parameter plane. From this figure, we cansee that the parameter r should be sufficiently large to show oscillatory behavior. Moreover, numericallywe have checked the direction and stability of Hopf-bifurcation and found that the Hopf-bifurcation isonly supercritical. We use Matlab 2017b software to produce all the figures ((4) to (8)) in this section.14 f Prey (a) Bifurcation diagram of prey w.r.t. f . f Predator (b) Bifurcation diagram of predator w.r.t. f .Figure 5: Bifurcation diagrams for Model (3.1) with respect to f . Here, f varies from 0 to 0 .
75, with the initial condition[ x (0) , y (0)] = [0 . , . r = 1 . c = 0 . d = 0 . d = 0 . a = 0 . α = 0 .
5, and m = 0 . r and f are r c = 0 . f c = 0 . c Prey (a) Bifurcation diagram of prey w.r.t. c . c Predator (b) Bifurcation diagram of predator w.r.t. c .Figure 6: Bifurcation diagrams for Model (3.1) with respect to c . Here, c varies from 0 to 35, with the initial condition[ x (0) , y (0)] = [0 . , . r = 0 . f = 1, d = 0 . d = 0 . a = 0 . α = 0 .
5, and m = 0 . r = d maα = 0 . c = c r = 5 . Prey
Predator E * = (0.2, 3.4) (a) Stable limit cycle oscillation occur around the co-existence equilibrium E ∗ , for f = 1. Prey
Predator E * = (0.2, 2.8) (b) Coexistence equilibrium is stable for f = 1 . r = 2, c = 0 . d = 0 . d = 0 . a = 0 . α = 0 .
5, and m = 0 .
01, corresponding to the Model (3.1).Figure 8: Hopf-bifurcation curve in r − f parameter plane for the Model (3.1). Other parameter values are fixed at c = 0 . d = 0 . d = 0 . a = 0 . α = 0 .
5, and m = 0 . . Discussion One of the central topics in ecology and evolutionary biology is to understand the variety of mech-anisms which influence the fitness and survivability of the population (Pough, 1989). In the literature,most of the studies only concentrated on the lethal effects of predator in prey-predator interaction. How-ever, some recent studies showed that apart from the lethal effects, there are some non-lethal effects,whose impact is equally important as of the previous one (Preisser and Bolnick, 2008). Among such ef-fects, fear of predation has an important role on population fitness and survivability in the prey-predatorsystem. Even such non-lethal effects are not restricted to affect in a single generation or in a particularseason, it can be carried-over over generations or within generation (O’Connor et al., 2014). Therefore,in the present study, we consider both the cost of predation fear and its carry-over effects in populationmodel.First, we developed and analyzed a single species (e.g., prey) population model, by incorporatingnon-lethal effects of predator by considering constant predator population and without explicit preda-tor dynamics. We incorporate the non-lethal effects in form of birth rate reduction due to the fear ofpredation. Moreover, we consider that such non-lethal effects can be carried within or over generations.We provide detailed analysis of our single species model, by neglecting direct predation, however, thequalitative properties of the model will be the same, if we consider predation followed by the Hollingtype I functional response. We derive the global dynamical properties of our proposed single speciesModel (2.2). For the single species Model (2.2), our main study objective is to find the different growthdynamics due to the cost of predation fear and its carry-over effects. More specifically, our goal is toinvestigate the occurrence of Allee effects due to such non-lethal effects of predator. From our anal-ysis, we can see that such non-lethal effects can be a cause of generating Allee effects and our modelshows three types of growth dynamics; namely, weak Allee dynamics, strong Allee dynamics and logisticdynamics, depending on the restrictions of model parameters, which is summarized in the Table (1).The above results not only relate the Allee effects and predation fear mechanism, but also answer thefirst objective raised in the introduction section. Moreover, our system shows both saddle-node andtranscritical bifurcations depending on the parameter values.Next, we include the explicit dynamics of predator population, where predation follows Holling type I (Holling, 1959) functional response. We derive the basic dynamical properties of our proposed model,existence and local stability conditions of each equilibria. Unlike the previous studies (Wang et al.,2016; Sasmal, 2018), our mathematical and numerical results show that the cost of predation fear andits carry-over effect, affect the prey-predator interactions in many ways, even if predation is followed bythe Holling type I functional response. If the birth rate of prey is small enough, then non-lethal effectsparameters have no effects on the system stability at coexistence equilibrium. However, if the birth rateof prey is large enough to support oscillation, then non-lethal effects parameters have great impact onthe system stability. System can be stable if non-lethal effects parameters are less than some thresholddensity, which is discussed in Theorems (3.2) and (3.3). Moreover, at the coexistence equilibrium, equi-librium density of prey does not depend on any of the non-lethal effects parameters, however, predatordensity decreases as we increase the parameter associated with cost of predation fear ( f ). On the otherhand, the effect of carry-over parameter ( c ) is opposite, i.e., predator density increases as we increasethe parameter ( c ), when population are at the coexistence steady state. The above discussion fulfillsthe third objective given in the introduction section regarding the role of non-lethal effects on systemstability at the coexistence equilibrium. 17e provide the global stability conditions of each equilibria in Theorems (3.4) and (3.5). As weincrease the birth rate parameter ( r ), our Model (3.1) experiences two bifurcations of equilibrium anda Hopf-bifurcation at positive equilibria. Existence of Hopf-bifurcation, its direction and stability atthe interior equilibrium are discussed in Theorems (3.6) and (3.7). From our analysis, we can see thatthe unique coexistence equilibrium is globally asymptotically stable if the birth rate of prey is not largeenough to support the oscillation. One of the interesting result is that our model system supports oscilla-tion, even for the Holling type I functional response, which can not be observed in classical prey-predatormodel with Holling type I functional response. Even oscillation cannot be observed in the presence ofonly growth rate reduction due to predation fear and with type I functional response (Wang et al.,2016; Sasmal, 2018). Therefore, the phenomenon ’paradox of enrichment’ (Rosenzweig, 1971; Gilpin andRosenzweig, 1972) can be observed in our model. However, both analytical and numerical results suggestthat such phenomenon can be ruled out by choosing suitable values of cost of fear and/or carry-overeffects. This answers the second and fourth questions listed in the introduction section regarding theglobal dynamical behavior and the phenomenon ’paradox of enrichment.’In the present study, we split the logistic dynamics into birth - death to incorporate the cost of fearonly in the birth rate. However, some theoretical study showed that due to the complexity of ecosystem,cost of fear can affect in many ways, like, it may increase the adult death rate, intra-specific competition,etc. (Zanette et al., 2011; Clinchy et al., 2013; Cresswell, 2011). Therefore, one may consider the fearand its carry-over effect directly to the logistic growth dynamics, when some experimental evidencesare available in future. We consider the simplest type I functional response in the presence of predatorpopulation. However, a more complicated functional response, like, Holling type II or III functionalresponse can be a mechanism for generating Allee effect in prey due to their predation satiation properties(Gascoigne and Lipccius, 2004). Moreover, two or more Allee effects can occur simultaneously in thesame population, and this is known as double or multiple Allee effects (Berec et al., 2007). Therefore,it may be interesting to see how multiple Allee effects occur due to fear and its carry-over effects whenpredator follows Holling type II or III functional response. Furthermore, it may be interesting to studyhow different dynamics can be observed by considering fear and its carry-over effects for other models,for example, Lotka-Volterra model, Beverton-Holt model etc.
6. ProofsProof of the Theorem (2.1)
Proof.
The local stability of x , x and x for Model (2.2) can be determined by the sign ofΥ( x i ) = Θ( x i ) + x Θ (cid:48) ( x i ) , where Θ (cid:48) ( x i ) = rcfy (1+ cx i + fy ) − d . Now, Υ( x ) = Θ(0) = r fy − d . Thus, the equilibrium x is locallyasymptotically stable if r < d (1 + f y ).For the first part of the Theorem (2.1), Φ( x ) > , ∀ x >
0, and consequently, Θ( x ) < , ∀ x >
0. As(2.2) is a scalar differential equation, x attracts every solution, i.e., x is globally asymptotically stable.Moreover, Θ( x ) can be written asΘ( x ) = − Φ( x )1 + cx + f y = ( x − x )( x − x )1 + cx + f y .
18t is easy to obtain, Θ (cid:48) ( x ) = − ( x − x )1+ cx + fy ( <
0) andΘ (cid:48) ( x ) = ( x − x )1+ cx + fy ( > . for x (cid:54) = x (as x > x for the case when Φ( x ) = 0 has real roots).Thus, Υ( x ) < x ) >
0. Therefore, the equilibrium x is always locally asymptoticallystable when it exists, whereas the equilibrium x is always unstable.For part 2( i ), x is unstable as r > d (1 + f y ). Moreover, Θ( x ) < ∀ x > x and Θ( x ) > ∀ x < x . As equation (2.2) is a scalar differential equation, then x attracts every solution in R + , i.e., x is globally asymptotically stable.For part 2( ii ), also the extinction equilibrium x is unstable and the positive equilibrium x is locallyasymptotically stable (as Θ (cid:48) ( x ) < x ) < ∀ x > x , andΘ( x ) > ∀ x < x , and therefore x attracts every solution in R + .For part 2( iii ), since ∆ = 0, then Ω must be positive, i.e., r < d (1 + f y ), which is the local stabilitycondition of x . Here, Φ( x ) = 0 has an identical real positive root, and Θ (cid:48) ( x ) = 0. Therefore, we can’tapply eigenvalue approach for the local stability. However, Φ( x ) > x > x (cid:54) = x . Therefore,Θ( x ) < ∀ x > x (cid:54) = x . Therefore, x is locally asymptotically stable as it attracts every solutionin R + with x < x and x is a saddle as it attracts every solution, starting at x > x , but repels x < x .For the final part, both the equilibria x and x are locally asymptotically stable. Here, Θ( x ) < ∀ x > x and 0 < x < x , Θ( x ) > ∀ x < x < x . Thus the basin of attraction of the extinctionequilibrium x is [0 , x ), and for the positive equilibrium x is R + \ [0 , x ). Proof of the Lemma (3.1)
Proof. As dxdt (cid:12)(cid:12)(cid:12) x =0 = 0 and dydt (cid:12)(cid:12)(cid:12) y =0 = 0 for any x ≥ y ≥
0, then x = 0 and y = 0 are invariantmanifolds, respectively. Due to the uniqueness of solution the set R is positively invariant for the Model(3.1). Moreover, dxdt ≤ rx (1 + cx )1 + cx + f y − d x − d x ≤ rx − d x − d x . By the comparison theory we can prove thatlim t →∞ sup x ( t ) ≤ r − d d . Define w ( t ) = x ( t ) + α y ( y ), then dwdt = rx (1+ cx )1+ cx + fy − d x − d x − mα y< rx − d x − d x − m ( w − x )= ( r − d + m ) x − d x − mw. Then similar to the proof of the theorem (2 .
1) in (Sasmal and Takeuchi, 2020), our results follow.19 roof of the Theorem (3.1)
Proof.
The characteristic equation at the interior equilibrium E ∗ is given by λ − Ψ λ + Ψ Ψ = 0 , (6.1)whose roots will be real negative or complex conjugate with negative real parts if Ψ <
0, whereΨ = x ∗ (cid:104) rcfy ∗ (1+ cx ∗ + fy ∗ ) − d (cid:105) Ψ = x ∗ (cid:104) rf (1+ cx ∗ )(1+ cx ∗ + fy ∗ ) + a (cid:105) ( > = aαy ∗ ( > . Now, Ψ < rcf y ∗ (1 + cx ∗ + f y ∗ ) < d , i.e., iff r < d ( aα + cm + aαf y ∗ ) a α cf y ∗ . Therefore, E ∗ is locally asymptotically stable if r < d ( aα + cm + aαfy ∗ ) a α cfy ∗ . The local stability of othertwo equilibria are similar, and hence not discussed here.Moreover, rcfy ∗ (1+ cx ∗ + fy ∗ ) − d < , ⇔ d a α f y ∗ + (cid:2) d aαf ( aα + cm ) − ra α cf (cid:3) y ∗ + d ( aα + cm ) > , which is always true if (for any real positive y ∗ ) (cid:2) d aαf ( aα + cm ) − ra α cf (cid:3) − d a α f ( aα + cm ) < , ⇒ r < d ( aα + cm ) aαc . Proof of the Theorem (3.2)
Proof.
We have, rcfy ∗ (1+ cx ∗ + fy ∗ ) − d < rcf y ∗ − d and rcf y ∗ − d < , ⇔ f rc (cid:2) aαr c ( aα + cm ) − rc ( aα + cm )( d aα + d m ) − aαd ( d aα + d m ) (cid:3) < a αd [ aαd + rc ( aα + cm )] . (6.2)Thus, the condition (6.2) is a sufficient condition for the local stability of E ∗ .Left hand side of the inequality (6.2) will be positive if r > c ( aα + cm )( d aα + d m )+ √ c ( aα + cm )( d aα + d m )[ c ( aα + cm )( d aα + d m )+4 a α d ]2 aαc ( aα + cm ) = r c . (6.3)When r > r c , then the sufficient condition for the local stability of E ∗ is f < a αd [ aαd + rc ( aα + cm )] rc [ aαr c ( aα + cm ) − rc ( aα + cm )( d aα + d m ) − aαd ( d aα + d m )] = f c . (6.4)20 roof of the Theorem (3.3) Proof.
We have, rcfy ∗ (1+ cx ∗ + fy ∗ ) − d < , ⇔ d a α f y ∗ + (cid:2) d aαf ( aα + cm ) − ra α cf (cid:3) y ∗ + d ( aα + cm ) > , a sufficient condition, that the above inequality holds true for any positive real y ∗ , is (cid:2) d aαf ( aα + cm ) − ra α cf (cid:3) − d a α f ( aα + cm ) < , ⇔ c ( raα − d m ) < d aα. (6.5)Thus, the condition r < d maα is a sufficient condition for the local stability of E ∗ .Since the left hand side of the inequality (6.5) is positive if r > d maα , and the sufficient condition forthe local stability of E ∗ is c < d aαraα − d m = c r . (6.6) Proof of the Theorem (3.4)
Proof.
The global stability of E in r ∈ (0 , d ) follows from lemma (3.1) and Theorem (3.1). Moreover,when r ∈ (cid:0) d , d + d maα (cid:1) , there exists only two equilibria E and E in R , and hence there can not be anyperiodic orbit in R , which implies that every solution will converge to any one of E and E . However,when r ∈ (cid:0) d , d + d maα (cid:1) then E is a saddle (repelling) and every solution will approach to E . Therefore,from local stability condition, E exists and is globally asymptotically stable if r ∈ (cid:0) d , d + d maα (cid:1) . Proof of the Theorem (3.5)
Proof.
A sufficient condition for the local stability of the positive equilibrium is r < d ( aα + cm ) aαc . Since, d c < d ( aα + cm ) aαc , from the previous Theorems (3.1) and (3.4), to show the global stability of E ∗ , it issufficient to prove that for r ∈ (cid:0) d + d maα , d c (cid:1) , there is no periodic orbit in { ( x, y ) | x > , y > } .By taking the Dulac function D ( x, y ) = xy for the System (3.1), we havediv (cid:12)(cid:12)(cid:12) ( D dxdt ,D dydt ) = ∂∂x (cid:0) D ( x, y ) dxdt ( x, y ) (cid:1) + ∂∂y (cid:16) D ( x, y ) dydt ( x, y ) (cid:17) = y (cid:104) rcfy (1+ cx + fy ) − d (cid:105) = y (1+ cx + fy ) (cid:2) − d f y + { rcf − d f (1 + cx ) } y − d (1 + cx ) (cid:3) . It is easy to see that [ rcf − d f (1 + cx )] − d f (1 + cx ) < r < d c .Therefore, div (cid:12)(cid:12)(cid:12) ( D dxdt ,D dydt ) < { ( x, y ) | x > , y > } if r < d c . Then by the Dulac-Bendixsontheorem (Perko, 2013), there is no periodic orbit in { ( x, y ) | x > , y > } for System (3.1) if r < d c .Moreover, E ∗ is the only stable equilibrium in { ( x, y ) | x > , y > } if d + d maα < r < d c . Hence everypositive solution will tend to E ∗ . 21 roof of the Theorem (3.7) Proof.
To find the stability and direction of Hopf-bifurcation, we calculate the 1st Lyapunov coefficient.Let u = x − x ∗ and v = y − y ∗ , then the System (3.1) becomes dudt = r ( u + x ∗ )(1+ c ( u + x ∗ ))1+ c ( u + x ∗ )+ f ( v + y ∗ ) − d ( u + x ∗ ) − d ( u + x ∗ ) − a ( u + x ∗ )( v + y ∗ ) := g ( u, v ) dvdt = aα ( u + x ∗ )( v + y ∗ ) − m ( v + y ∗ ) := h ( u, v ) . Now, considering the Taylor’s series expansion at ( u, v ) = (0 ,
0) up to 3rd order, we have dudt = g u u + g v v + g ( u, v ) , dvdt = h u u + h v v + h ( u, v ) , (6.7) g ( u, v ) and h ( u, v ) are the higher order terms of u and v , given by g ( u, v ) = g uu u + g uv uv + g vv v + g uuu u + g uuv u v + g uvv uv + g vvv v ,h ( u, v ) = h uu u + h uv uv + h vv v + h uuu u + h uuv u v + h uvv uv + h vvv v , where g u = − d x ∗ + rcfx ∗ y ∗ (1+ cx ∗ + fy ∗ ) , g v = − (cid:104) ax ∗ + rfx ∗ (1+ cx ∗ )(1+ cx ∗ + fy ∗ ) (cid:105) , g uu = − d − rcfy ∗ ( − cx ∗ − fy ∗ )(1+ cx ∗ + fy ∗ ) ,g uv = rcfx ∗ (1+ cx ∗ − fy ∗ )(1+ cx ∗ + fy ∗ ) , g vv = − rx ∗ (1+ cx ∗ ) f (1+ cx ∗ + fy ∗ ) , g uuu = rc fy ∗ ( − cx ∗ − fy ∗ )(1+ cx ∗ + fy ∗ ) ,g uuv = − rcf ( − c x ∗ + f y ∗ − cx ∗ fy ∗ )(1+ cx ∗ + fy ∗ ) , g uvv = − rcf x ∗ (2+2 cx ∗ − fy ∗ )(1+ cx ∗ + fy ∗ ) , g vvv = rx ∗ (1+ cx ∗ ) f (1+ cx ∗ + fy ∗ ) , and h u = aαy ∗ , h v = 0 , h uu = 0 , h uv = aα, h vv = 0 , h uuu = 0 , h uuv = 0 , h uvv = 0 , h vvv = 0 . Here all the partial derivatives are calculated at the bifurcation point, i.e., ( u, v ) = (0 , U = (cid:20) g u g v h u h v (cid:21) U + F ( u ) , where U = ( u, v ) T and F = ( f ( u, v ) , g ( u, v )) T = (cid:0) g uu u + g uv uv + g vv v + g uuu u + g uuv u v + g uvv uv + g vvv v , h uv uv (cid:1) T . Now, Hopf-bifurcation occurs when g u = 0, i.e., at the Hopf-bifurcation point, the eigenvalue will bepurely imaginary, which is given by iω , where ω = √− g v h u . Eigenvector corresponding to this eigenvalue iω is given by ¯ v = ( g v , iω ) T . Now, we define Q = ( Re (¯ v ) , − Im (¯ v )) = (cid:20) g v iω (cid:21) . Now, let U = QZ or Z = Q − U , where Z = ( z , z ) T . Therefore, under this transformation, the system is reduced to˙ Z = (cid:18) Q − (cid:20) g u g v h u h v (cid:21) Q (cid:19) Z + Q − F ( QZ ) . This can be written as (cid:20) ˙ z ˙ z (cid:21) = (cid:20) − ωω (cid:21) (cid:20) z z (cid:21) + (cid:20) F ( z , z ) F ( z , z ) (cid:21) , where F ( z , z ) and F ( z , z ) are given by 22 ( z , z ) = g v (cid:2) g uu g v z − ωg uv g v z z − g v h u g vv z + g uuu g v z − ωg uuv g v z z − g v h u g uvv g v z z − g v h u ωg vvv z (cid:3) F ( z , z ) = g v h uv z z . The direction of Hopf-bifurcation is determined by the sign of the 1st Lyapunov coefficient, which isgiven by L := (cid:104) ∂ F ∂z + ∂ F ∂z ∂ z + ∂ F ∂ z ∂z + ∂ F ∂ z (cid:105) + ω (cid:104) ∂ F ∂z ∂z (cid:16) ∂ F ∂ z + ∂ F ∂ z (cid:17) − ∂ F ∂z ∂z (cid:16) ∂ F ∂ z + ∂ F ∂ z (cid:17) − ∂ F ∂ z ∂ F ∂ z + ∂ F ∂ z ∂ F ∂ z (cid:105) . We use the maple software to simplify the expression of L , which is given as follows: L := 3 g v g uuu [(1 − g uv z ) + 2 z ( g uvv z w − g uuv z g v )] + 3 wg vvv z [2 w ( h u g uvv z − wg uuv z ) + g uv h u ]+ g uvv w [ w (1 − g uv z ) + 2 z ( g v g uu − h u g vv ) + 2 z w ( wg uvv z − g v z g uuv )]+ 2 g uuv z (cid:2) g v ( wg uuv z − g uu ) − w ( g vv + g v z g uvv ) (cid:3) + g uv ( h u g vv + wg v g uuv z − g v g uu ) . Now by (Perko, 2013), Hopf-bifurcation is supercritical if
L <
L > Acknowledgement
YT’s research is supported by Aoyama Gakuin University research grant “Ongoing Research Support”and Japan Society for the Promotion of Science “Grand-in-Aid 20K03755.”