A. K. Misra, S. Pal & Y. Kang
Species extinction rate is unprecedentedly growing over time. According to the International Union of Conservation of Nature (IUCN), 160 species have gone extinct in the last decade [1], which includes mammals, reptiles, birds and plants. Furthermore, 18% of extant vertebrates have been declared vulnerable [2]. Atwood et al. [7] concluded from their study that herbivores perceive elevated predation risk among mammals, birds and reptiles. Some theoretical studies with the help of mathematical modeling have been done to preserve endangered species. Oliveira and Hilker [32] explored bio-control approach by introducing disease in invasive predator species. Disease weakens the predators and the augmented death results in limiting the predation, which in turn will allow the endangered species to recover. Numfor et al. [31] studied optimal control strategy of trapping and culling of invasive predators along with bio-control strategy to conserve endangered species. As every organism has fair contribution and specific role in shaping and preserving the healthy ecosystems in which they dwell, continuous species loss has influenced the ecology of our planet profoundly. All the species in terrestrial or aquatic ecosystem occupy some trophic level in the food chain. Often, extinction of lower trophic species leads to secondary extinction. Predation or feeding serves as the bridge between trophic levels, which provides the pathway for energy to flow from lower trophic levels to higher trophic levels. The entire idea of energy flow mechanism is based on the bottom-up approach, i.e. the cumulative resources, like food and habitat available for lower trophic level ultimately determine the fate of those species occupying higher trophic levels. This concept of bottom-up approach has been considered as one of the main tenets of ecology. However, to answer the famous question "why is the Earth green?", ecologists agreed that influence of top-down approach along with bottom-up approach is also pervasive [20,38]. In this top-down approach, the role of predators becomes significant in shaping the ecosystem. Carnivores weeding out weak, slow and dying animals, which belong to their prey community, and thus keep prey population in check. The presence of predators restrains herbivores from accessing and over-utilizing the plant resources they feed upon. Thus, predation is not only beneficial for predators, but also keeps the entire ecosystem healthy. Therefore, ecologists come out with the concept of introducing top predator as one of the most plausible techniques to improve ecological functioning. Sergio et al. [44] reviewed the role of top predators in bio-diversity restoration in reference of various ecosystems. Baker et al. [8] proposed an ensemble modeling method to study the potential outcomes of keystone predator reintroduction in an ecosystem.
An experiment of repatriation of wolf in Yellowstone national park in the USA comes out with resounding success in reviving the degrading ecosystem. In 1995, eight wolves were reintroduced in Yellowstone national park with the expectation of restoring the continuously deteriorating ecosystem. Reintroduction of wolves brought back lots of ecological benefits for the whole national park ecosystem [39]. As predation by wolves limits the population of deer and elk, some endangered 2450068-2 plant species were observed in more patches, colony of beaver species increased, more songbirds were being observed as canopy increased, more scavengers were seen as carcasses increased. Some barren land came to life as grazing decreased by some amount. Besides that, more interestingly, it is observed that the behavior and grazing pattern of elk have changed in the presence of wolves. Very recently, in 2018, the reintroduction of wolves in Isle Royale Island also got reverberating success in bringing back the degrading ecosystem to its healthy state [48].
Mathematical modeling provides a platform to understand the dynamic process involved in ecology and are often useful to make practical predictions and derive insightful conclusions. Starting from the seminal work of Lotka-Volterra [30,54], theoretical study of predator-prey interaction and food web system traverse a long way with the help of mathematical modeling and has proven to be a persuasive alternative for time consuming and often risky field experiments. Development of mathematical models with the incorporation of species specific traits engaged many researchers working in the field of ecological modeling. Allee, in 1930s, put forward his observation that in many species, the growth and population density are positively correlated. Biological phenomena like difficulty in mating and reduced anti-predator defense lead to the Allee effect. Numerous predator-prey systems have been studied considering strong Allee effect [3,15,55,56,62,63] for its ecological significance. Along with the modified growth term, predator functional response plays a key role in modulating predator-prey dynamics. Predation usually involves searching of food entity and food sharing. Ratio-dependent functional response, that stands on the ratio of prey and predator population instead of depending on only prey population, better captures the mutual interference among predators [13,19]. Interestingly, for ratio-dependent predator-prey system, both the populations may extinct even in presence of stable co-existence equilibrium in the system. Many researchers studied the predator-prey dynamics considering functional response to be ratio-dependent [15,16,21,23,53,60]. Besides the consumptive effect of predator, the growth of prey species is also influenced by the non-consumptive effect. That is, only the presence of predator brings substantial psychological and behavioral changes in prey species [14,46,49,59]. The change in foraging behavior and spending more effort and time in vigilance to counter the fear of predation attributes to reduction in the growth rate of species. A significant number of modeling-based studies have been done to see the impact of fear on the dynamics of different predator-prey system [18, 28, 40-42, 51, 57]. Cong et al. [12] and Verma et al. [52] studied the role of fear in a three-species food chain model. Panday et al. [36] observed that fear induces trophic cascading in a three species food chain model.
An ecosystem is seldom stable, rather it is a very dynamic system and is always in recovering phase. However, sometimes, a particular entity of an ecosystem, i.e. a species becomes badly affected mainly due to over-exploitation and environmental factors. The affected species may even reach at the verge of extinction. In this paper, we provide a theoretical study to avoid the extinction of a species that occupy a 2450068-3 relatively lower trophic level in a food chain. We conceptualize the introduction of apex predator as a species conservation tool.
Specific research questions, we intend to address in this study, are as follows:
• If the density of a species depleted below a certain threshold, i.e. viable population size, which triggers the chance of extinction of the species, then how the extinction can be avoided and density of the species can be restored naturally? • In modeling phenomena, we incorporate the impact of predation induced fear along with direct killing. How these fear parameters mediate the dynamics of the proposed systems?
Suppose x is the density of our target species, which takes the place of prey in a predator-prey interaction. y is the density of predator (mesopredator) population, which totally depends on the considered prey species for food. In previous studies, predator functional response, i.e. per capita predator's food consumption per unit time, was thought to be a function of prey density only (Holling type I, II, and III functional responses), which covers a huge volume of literature. However, over the course of time, with ecological justifications, it is believed that predator density also has a role to play in the predator's functional response. Thereafter, ratiodependent functional response, a particular form of predator density-dependent response, where food consumption by a predator per unit time is a function of ratio of prey density to predator abundance, becomes popular among many researchers.
Outcomes of many field and laboratory experiments support the consideration of ratio-dependent functional response [5,6,9]. Kuang and Baretta [27] first analyzed global qualitative behavior of a ratio-dependent predator-prey system in a systematic way. Thereafter, Hsu et al. [23] contributed with a more detailed study of global qualitative behavior and answered many open questions left unanswered by Kuang and Baretta. Xiao and Ruan [58] also put light on global behavior in the neighborhood of origin of ratio-dependent predator-prey system. They have investigated different kinds of topological structures in the neighborhood of origin. In many terrestrial and aquatic ecosystems, the growth rate of some species is observed to follow Allee dynamics. At smaller species density, i.e. below a certain threshold population level, net population growth becomes negative. Few studies have been carried out by considering ratio-dependent functional response along with the incorporation of Allee effect in the growth rate of species [3,10,26,43]. In particular, Sen et al. [43] explicitly incorporated the Allee factor along with logistic growth term of prey in a predator-prey system with ratio-dependent functional response. In this study, they compare the dynamics of ratio-dependent predator-prey model with and without Allee effect. They concluded that chance of persistent oscillation of species density ceases with the consideration of Allee effect in the ratio-dependent predator-prey system. Surprisingly, the study of non-consumptive effect of predator over prey density (effect of predation fear), where prey's growth rate subjected 2450068-4
to strong Allee effect along with ratio-dependent predator's functional response, is overlooked. Here, first we consider a two species system, where we represent the interaction of two species by a classical predator-prey model with the ratiodependent Michaelis-Menten-type functional response [24,25]. Furthermore, the non-consumptive impact of predators, namely, the induction of fear, elicits significant behavioral and physiological alterations in prey species. Predation-induced fear not only contributes to diminished reproductive output [59], but also affects foraging behavior and adult survival [4,11,47]. Consequently, the fear of predation influences the inherent reproductive capacity of prey populations. Thus, the reproductive rate is modulated by a declining function of predator abundance and the intensity of predation-induced fear [35,50,61]. Let h(k 1 , y) represent this decreasing function, where k 1 denotes the intensity of fear and y signifies predator abundance.
The biologically pertinent assumptions underlying the formulation of such a function are elucidated as follows:
The simplest function with these properties is h(k 1 , y) = 1 1+k1y , which is also considered by Wang et al. [57] and many authors [12,28,33,34] to include the impact of predation induced fear in the growth of prey species. With the incorporation of fear factor and Allee effect in the growth equations of prey, the predator-prey system takes the following form:
where x(0) > 0 and y(0) > 0. As growth equations of both prey and predator are undefined at (x, y) = (0, 0), we redefine growth rates dx dt = 0, dy dt = 0, at (x, y) = (0, 0). Such modification was first proposed by Xiao and Ruan [58], and thereafter considered in many studies.
It can be easily shown that for 0 < x(0) < A, lim t→∞ x(t) = 0, which in turn implies y(t) → 0 as t → ∞. Therefore, for our proposed model system, we consider x(0) > A. To understand the relationship between population size of a species and its chances of extinction, ecologists came up with the idea of viable population size [22,45]. It is the minimum population size range at which preservation method can be successfully accomplished. Prior to applying any conservation strategy on any species, it is very important to predict its viable population size. Let x c be that critical limit of population size, below which the species is in danger of extinction. This critical value varies from species to species and also depends on various factors such as nature of ecosystem, habitat and environmental factors.
When the population of x-species is large enough, i.e. x > x c , where x c is an arbitrary chosen threshold population size of x-species, the dynamics is governed by system (2.1). If the population size of x-species depletes and enters into the interval A < x ≤ x c , where A is the Allee threshold, then to protect the x-species from further diminution, we consider that top predator species that feeds upon yspecies (mesopredator) only but not on x-species, is being introduced in the habitat. With the introduction of top predator, the dynamics of 3-species system is assumed to be governed by the following system of differential equations:
(2.2)
Initial conditions are x(0) > 0, y(0) > 0, and z(0) > 0. Again, the growth equations are defined to be dx dt = 0, dy dt = 0, and dz dt = 0, at (x, y, z) = (0, 0, 0). Here, along with the impact of direct predation of top predator, the non-consumptive effect is also taken into account in governing the dynamics of mesopredator and hence prey species. In addition to impeding growth, fear also impacts the effective predation by mesopredators. A field investigation by Gordon et al. [17] furnished compelling evidence of the suppression of foraging and predatory behavior among mesopredators (e.g. feral cats) in the presence of apex predators (e.g. dingoes), thereby mitigating the perceived predation risk encountered by small prey, such as desert rodents. Moreover, within the Yellowstone National Park ecosystem, the reintroduction of apex predators like wolves has brought about substantial shifts in the behavior and grazing patterns of elk [39]. Consequently, in our proposed model, both the growth and predation terms concerning mesopredators are subject to multiplication by a factor of (1/1 + k 2 z), representing a diminishing function of the fear parameter (k 2 ) and the density of top predators (z). Description of all the parameters that appeared in system (2.2) is mentioned in Table 1. It is to be noted that, we define x-species as prey, y-species as mesopredator and z-species as top predator, throughout the paper. We intend to study the effect of introduction of top predator on the density of our targeted prey species and explore the whole dynamics exhibited by systems (2.1) and (2.2). Now, we show that solution of system (2.1) satisfies the positivity and boundedness criteria. From system (2.1), we have
-δ 1 du .
How can we avoid the extinction of any species naturally? Density-dependent death rate of top predator Therefore, x(0) > 0 and y(0) > 0 imply x(t) ≥ 0, y(t) ≥ 0, ∀ t ≥ 0. Hence, the solution trajectories initiated from positive quadrant of x -y plane stay in the positive quadrant.
To show boundedness of the system, consider
where
. Now, we discuss two cases.
Case I. x(0) ∈ (0, K). We intend to show x(t) ≤ K, ∀ t ≥ 0. On the contrary, suppose there exist two positive values of t; T 1 and T 2 such that x(T 1 ) = K and x(t) > K, ∀ t ∈ (T 1 , T 2 ). Then for all t ∈ (T 1 , T 2 ),
which is contradictory to our assumption. Therefore, x(t) ≤ K for all t > 0. Case II. x(0) > K.
As F (x(t), y(t)) ≤ 0 for x(t) ≥ K, therefore whenever x(t) ≥ K,
Combining two cases, we get
Again, system (2.1) implies
Using Gronwall's inequality, we have
).
As δ 1 > 0, so for sufficiently large time, we get the inequality
takes any positive value. Therefore, y(t) is also bounded. In a similar manner, for system (2.2), we can have
Therefore, with similar argument as above, x(t) and y(t) are bounded as for large time, (x(t)+ 1 θ1 y(t)) ≤ γ δ1 + , takes any positive value. Now, dz dt ≤ θ 2 α 2 y max z-δ 2 z 2 implies z(t) ≤ max{z(0), θ2α2ymax δ2 }, y max is the upper bound of y.
3. Dynamics of the System (2.1)
Along with trivial equilibrium E 0 (0, 0), the predator-free equilibria E 1 (A, 0) and E 2 (K, 0) are always feasible. Solution of the following two algebraic equations gives feasible coexistence equilibrium E * (x * , y * ):
From the second equation, we have
Using this value of y in first equation, we get
Simplifying, we get
where
Therefore, system (2.1) has two co-existence equilibria if θ 1 α 1 -δ 1 > 0, a 1 < 0 and a 2 1 > 4a 2 . These two equilibria collide when a 2 1 = 4a 2 and vanish for a 2 1 < 4a 2 . These conditions raise the possibility for system (2.1) to exhibit saddle-node bifurcation, provided θ 1 α 1 -δ 1 > 0 and a 1 < 0 hold.
Trivial equilibrium E 0 (0, 0) is non-hyperbolic attractor for all parameter values, which can be proved as per [15,Lemma 3]. The equilibrium E 1 (A, 0) is always unstable, whereas the stability of equilibrium E 2 (K, 0) depends on the sign of θ 1 α 1δ 1 . Equilibrium E 2 is stable only when the maximum per capita growth rate of mesopredator is negative, i.e. θ 1 α 1 -δ 1 < 0. In that case, only prey population survive and mesopredator population extinct. E 2 becomes unstable when the sign of θ 1 α 1 -δ 1 changes to positive. Thus, system (2.1) shows transcritical bifurcation between E 2 (K, 0) and coexistence equilibrium, provided co-existence equilibrium exist. The Jacobian matrix at E * (x * , y * ) is
2450068-9
where h
, and
. Therefore, we can conclude that E * is locally stable only when
Also, it is very clear that if h (x * ) is negative, then the equilibrium E * is always stable.
Remark 1. From Eq. (3.3), it is to be noted that the x-component of the stable
, where a 1 < 0 and a 2 1 -4a 2 > 0. Now,
rδ1 > 0 and da2 dα1 = K r > 0. Therefore, it is evident that equilibrium density of prey species always follow decreasing trend with increasing k 1 and α 1 , whenever the equilibrium is feasible.
Here, first we show the occurrence of transcritical bifurcation by using Sotomayor theorem [37]. As E 2 (K, 0) changes its stability at θ 1 α 1 -δ 1 = 0. Therefore, we will look for transcritical bifurcation around E 2 (K, 0) with respect to the parameter α 1 . At E 2 , the Jacobian matrix of model system (2.1) is
At α 1 (= α * 1 ) = δ1 θ1 , the matrix J E2 has simple zero eigenvalue. The eigenvectors of J E2 and J T E2 corresponding to zero eigenvalue are V = (v 1 , v 2 ) T = ( -α1 r(K-A) , 1) T and W = (w 1 , w 2 ) T = (0, 1) T , respectively.
Consider, G = (g 1 , g 2 ) T , where
How can we avoid the extinction of any species naturally?
Now,
Therefore, according to Sotomayor theorem [37], the conditions for the existence of transcritical bifurcation are satisfied. Hence, transcritical bifurcation occurs at α 1 = α * 1 for system (2.1), that is as α 1 crosses α * 1 , the equilibrium E 2 changes its stability and emergence (or vanishing) of a stable interior equilibrium is observed.
As mentioned in Remark 1, two co-existence equilibria collide and vanish at (x, ỹ), where
. Corresponding x and ỹ are given by
and ỹ = θ1α1 δ1 x. As two equilibrium collide at (x, ỹ), discriminant of Eq. (3.3) becomes zero, which gives a critical value of k 1 , k1 say. Now, we check the conditions of Sotomayor theorem [37] for saddle-node bifurcation.
Let the eigenvectors corresponding to zero eigenvalue be V = (v 1 , v2 ) T and W = ( w1 , w2 ) T for J Ẽ and J T Ẽ , respectively. Then
Hence, by Sotomayor theorem, we affirm that system (2.1) undergoes saddle-node bifurcation around Ẽ = (x, ỹ) as k 1 crosses a critical value k1 .
Occurrence or termination of limit cycle with varying parameter is characterized by Hopf bifurcation. Choosing strength of fear k 1 as the bifurcation parameter, we 2450068-11 intend to show analytically that system (2.1) exhibits Hopf bifurcation at interior equilibrium E * (x * , y * ). Conditions for occurrence of Hopf bifurcation for system (2.1) are mentioned in the following theorem.
Suppose the eigenvalues of the Jacobian matrix for any neighboring point
, where
Transversality condition affirms the crossing of imaginary axis for the eigenvalues with nonzero speed. The transversality condition is satisfied if
Therefore, with these conditions, system (2.1) would undergo Hopf bifurcation at
4. Dynamics of System (2.2)
Positive solution of following three algebraic equations gives the co-existence equilibrium of system (2.2):
How can we avoid the extinction of any species naturally?
From (4.3), we have
Using this value of z in Eqs. (4.1) and (4.2), we get
)
The curve (4.6) approaches origin and always increases in the positive quadrant. On the other hand, the curve (4.5) passes through (A, 0) and (K, 0); increasing at (A, 0) and decreasing at (K, 0). Therefore, these two isoclines may have no intersection, 4.2. Stability and bifurcation analysis of system (2.2)
Now, the Jacobian matrix at co-existence equilibrium Ê = (x, ŷ, ẑ) is
where
The characteristic equation corresponding to the Jacobian matrix is
where According to the Routh-Hurwitz criterion, the roots of Eq. (4.8) lie on the left half of a complex plane if and only if B 1 > 0, B 3 > 0, and B 1 B 2 -B 3 > 0. Therefore, with these conditions, the equilibrium Ê = (x, ŷ, ẑ) is locally asymptotically stable.
Suppose at some critical value of α 2 (α * 2 say), the conditions B 1 > 0 and B 3 > 0 hold but B 1 B 2 -B 3 = 0. Then, the characteristic equation (4.8) can be written as
(4.9)
The above equation has two purely imaginary roots, say μ 1,2 = ±ι √ B 2 and a negative real root, say μ 3 = -B 1 .
How can we avoid the extinction of any species naturally?
Now, suppose at any point α 2 in the -neighborhood of α * 2 , μ 1,2 = γ 1 ± ιγ 2 . Putting this in the characteristic equation and separating real and imaginary parts, we get
As γ 2 = 0, from Eq. (4.11), we have
Using this value of γ 2 in Eq. (4.10), we get
Differentiating above equation with respect to α 2 and using the fact that γ 1 (α * 2 ) = 0, we get
Hence, transversality condition holds if
Therefore, with this condition, system (2.2) exhibits Hopf bifurcation at α 2 = α * 2 around interior equilibrium Ê. System (2.2) also shows saddle-node bifurcation, which can be proved following the same analysis carried out in Sec. 3.2.
Bogdanov-Takens bifurcation (BT-bifurcation) is a type of codimension-2 bifurcation. At BT-bifurcation point, the system has zero eigenvalue with multiplicity 2. Using the methods from Kuznetsov [29], we derive the conditions for BT-bifurcation at Ê of the model system (2.2).
To examine the conditions for BT-bifurcation of system (2.2) at Ê, first we use the transformations x = x + X, y = ŷ + Ȳ , and z = ẑ + Z in order to shift the origin at Ê = (x, ŷ, ẑ). Consequently, system (2.2) can be expressed in the subsequent form:
where
and 3 ,
How can we avoid the extinction of any species naturally?
Let P = [ Ũ1 Ũ2 Ũ3 ], then under the non-singular linear transformation
where the inverse of P is given by
Here, Ã20 = ṽ11 (l
1 200 ũ2 11 + l 1 020 ũ2 21 + l 1 002 ũ2 31 + l 1 110 ũ11 ũ21 + l 1 101 ũ11 ũ31 + l 1 011 ũ21 ũ31 ) + ṽ12 (l 2 200 ũ2 11 + l 2 020 ũ2 21 + l 2 002 ũ2 31 + l 2 110 ũ11 ũ21 + l 2 101 ũ11 ũ31 + l 2 011 ũ21 ũ31 ) + ṽ13 (-δ 2 ũ2 31 + θ 2 α 2 ũ21 ũ31 ), Ã11 = ṽ11 (2l 1 200 ũ11 ũ12 + 2l 1 020 ũ21 ũ22 + 2l 1 002 ũ31 ũ32 + l 1 110 (ũ 11 ũ22 + ũ12 ũ21 ) + 2l 1 101 (ũ 11 ũ32 + ũ12 ũ31 ) + 2l 1 011 (ũ 22 ũ31 + ũ21 ũ32 )) + ṽ12 (2l 2 200 ũ11 ũ12 + 2l 2 020 ũ21 ũ22 + 2l 2 002 ũ31 ũ32 + l 2 110 (ũ 11 ũ22 + ũ12 ũ21 ) + 2l 2 101 (ũ 11 ũ32 + ũ12 ũ31 ) + 2l 2 011 (ũ 22 ũ31 + ũ21 ũ32 )) + ṽ13 (-2δ 2 ũ31 ũ32 + θ 2 α 2 (ũ 21 ũ32 + ũ22 ũ31 )), Ã02 = ṽ11 (l 1 200 ũ2 12 + l 1 020 ũ2 22 + l 1 002 ũ2 32 + l 1 110 ũ12 ũ22 + l 1 101 ũ12 ũ32 + l 1 011 ũ22 ũ32 ) + ṽ12 (l 2 200 ũ2 12 + l 2 020 ũ2 22 + l 2 002 ũ2 32 + l 2 110 ũ12 ũ22 + l 2 101 ũ12 ũ32 + l 2 011 ũ22 ũ32 ) + ṽ13 (-δ 2 ũ2 32 + θ 2 α 2 ũ22 ũ32 ), B20 = ṽ21 (l 1 200 ũ2 11 + l 1 020 ũ2 21 + l 1 002 ũ2 31 + l 1 110 ũ11 ũ21 + l 1 101 ũ11 ũ31 + l 1 011 ũ21 ũ31 ) + ṽ22 (l 2 200 ũ2 11 + l 2 020 ũ2 21 + l 2 002 ũ2 31 + l 2 110 ũ11 ũ21 + l 2 101 ũ11 ũ31 + l 2 011 ũ21 ũ31 ) + ṽ23 (-δ 2 ũ2 31 + θ 2 α 2 ũ21 ũ31 ), 2450068-17 B11 = ṽ21 (2l 1 200 ũ11 ũ12 + 2l 1 020 ũ21 ũ22 + 2l 1 002 ũ31 ũ32 + l 1 110 (ũ 11 ũ22 + ũ12 ũ21 ) + 2l 1 101 (ũ 11 ũ32 + ũ12 ũ31 ) + 2l 1 011 (ũ 22 ũ31 + ũ21 ũ32 )) + ṽ22 (2l 2 200 ũ11 ũ1 + 2l 2 020 ũ21 ũ22 + 2l 2 002 ũ31 ũ32 + l 2 110 (ũ 11 ũ22 + ũ12 ũ21 ) + 2l 2 101 (ũ 11 ũ32 + ũ12 ũ31 ) + 2l 2 011 (ũ 22 ũ31 + ũ21 ũ32 )) + ṽ23 (-2δ 2 ũ31 ũ32 + θ 2 α 2 (ũ 21 ũ32 + ũ22 ũ31 )), B02 = ṽ21 (l 1 200 ũ2 12 + l 1 020 ũ2 22 + l 1 002 ũ2 32 + l 1 110 ũ12 ũ22 + l 1 101 ũ12 ũ32 + l 1 011 ũ22 ũ32 ) + ṽ22 (l 2 200 ũ2 12 + l 2 020 ũ2 22 + l 2 002 ũ2 32 + l 2 110 ũ12 ũ22 + l 2 101 ũ12 ũ32 + l 2 011 ũ22 ũ32 ) + ṽ23 (-δ 2 ũ2
32 + θ 2 α 2 ũ22 ũ32 ). Hence, in accordance with the center manifold theorem, center manifold exists for system (2.2), and it can be locally expressed in the following manner:
for and sufficiently small 1 and 2 . Thus, we have to calculate the center manifold for system (2.2). The system is restricted to the central manifold given as
Using the transformation
and rewriting system (4.14) in X1 and X2 , we get
where B 20 = B20 and B 11 = B11 + 2 Ã20 . The conditions derived for the occurrence of BT-bifurcation are summarized in Theorem 4.1.
Theorem 4.1. If B 20 and B 11 are nonzero, then system (2.2) displays a codimension 2 BT-bifurcation at the interior equilibrium Ê.
This section deals with numerical simulations to explore the underlying dynamics that systems (2.1) and (2.2) posses, with the help of MATLAB software and MAT-CONT package. Suppose initially the density of prey is large enough such that it
How can we avoid the extinction of any species naturally?
remains above the threshold value x c . Then the dynamics is governed by system (2.1) and the parameter values we consider for simulations are as follows:
First, we plot the equilibrium curve by varying fear parameter k 1 to capture the impact of fear (Fig. 2(a)). Supercritical Hopf bifurcation occurs at k 1 = 1.9956 (denoted by H) as first Lyapunov coefficient (l 1 ) is negative (l 1 = -1.828425 × 10 -7 ). Saddle-node bifurcation also is found to occur for system (2.1) at k 1 = 2.7194 (denoted by LP ). At LP , two branches of equilibrium curves collide and disappear. Continuation of limit cycles from the Hopf point has been plotted in Fig. 2(b). Stable limit cycles originate from Hopf point which further disappear at k 1 = 1.997, with the occurrence of Limit Point Curve (LPC). At LPC, saddle-node bifurcation of limit cycle occurs, i.e. one stable and one unstable limit cycle collide and vanish.
It is evident that prey population diminishes with the increment in fear parameter. Between Hopf point k 1 = 1.9956 and k 1 = 1.996, the dynamical variables show persistent periodic oscillations. Further, between k 1 = 1.996 and k 1 = 1.997, an unstable limit cycle surrounding stable limit cycle exists. Trajectories converge to stable limit cycle only when the initial point lies inside the outer unstable limit cycle. For initial points in the exterior of the unstable limit cycle, the trajectories collapse. If k 1 > 1.997, then species population goes to extinction. Therefore, the prey species, which shows more anti-predator behavior, is more vulnerable to extinction.
It is obvious that for system (2.1), whenever x(0) < A, solution trajectories converge to the trivial equilibrium E 0 (0, 0). However, this is not the necessary condition, i.e. even if x(0) > A, the solution trajectories may move toward the trivial equilibrium depending on the initial start. In Fig. 3, we have plotted the
How can we avoid the extinction of any species naturally? 1.5 2 2.5 3 3.5 4 4.5 5 0 0.5 1 1.5 2 2.5 3 Fig. 5. Bifurcation diagram in k 1α 1 parametric plane.
As we see that predation induced fear and rate of predation have the potential to modulate system's dynamics, two-parameter bifurcation diagram is plotted in k 1 -α 1 plane (Fig. 5). Black curve denotes transcritical curve, below which coexistence equilibrium is not feasible, only mesopredator-free equilibrium is stable. In the parametric region between transcritical and Hopf curve, one stable co-existence equilibrium exists, which loses its stability as Hopf-curve is crossed. In the above saddle-node curve, co-existence equilibrium loses its feasibility again. We can infer from the figure that for less fearful prey, both the species may sustain stably in the habitat for considerably large range of predation rate. However, if the impact of fear is large enough, then higher predation rate always trigger extinction of species. Therefore, for more fearful prey, extinction is more likely to occur. Also, from Fig. 5, it is noted that if the parameter values lie in the region between transcritical curve
and Hopf curve, then the threshold value x c determines whether the introduction of top predator should be initiated in the habitat or not. However, whenever the parameters lie above the Hopf curve, introduction of top predator is inevitable, whatsoever x c might be.
From the above discussion, we observe that, for higher strength of fear and magnitude of predation rate, prey species equilibrium density is reducing, and even the collapse of species population is also occurring. In such a situation, the idea of introduction of top predator which predates on mesopredator but not on the considered prey species is deployed. In Fig. 6, we first plot time series for system (2.1) up to t = 50 unit, taking k 1 = 1.8 and other parameter values are same as given in Eq. (5.1). We see that system stabilizes at (x * , y * ) = (632.531, 253.088). Then at t = 50 unit, we introduce top predator in the system, and we plot time series from t = 50 unit to t = 100 unit, for system (2.2). The parameter values we consider for the simulations of system (2.2) are given by
(
We observe that top predator introduction induces trophic cascading effect and system (2.2) stabilizes with very high density of prey species. In Fig. 7, we have plotted the phase portraits of systems (2.1) and (2.2) for k 1 = 1.8. It is evident that the basin of attraction of the co-existence equilibrium is broadening with the incorporation of the top predator. Suppose we choose the critical value x c = 400, then from Fig. 7(a), we can see that trajectories with initial start outside the basin of attraction of stable interior equilibrium, i.e. upper half of green curve, bound to cross x c and ultimately falls to origin. The instant it crosses x c , we will introduce top predator. With the introduction of top predator,
How can we avoid the extinction of any species naturally? the species extinction can be controlled if the solution trajectories lie inside the basin of attraction for the equilibrium of system (2.2). In Fig. 7(b), z(0) is taken to be z(0) = 5 for all the solution trajectories, i.e. we are considering that 5 unit of top predator is being introduced. It is also observed that the value of z(0) has not much influence on the basin of attraction. Even though for system (2.2) the basin of attraction of the stable co-existence equilibrium is being increased, the trajectories may also converge to origin if the initial point lie on the left of green curve (Fig. 7(b)).
Similarly, in Fig. 8, we plot time series for system (2.1) in time interval [0, 50] and for system (2.2) in [50,100], taking k 1 = 1.9965 and other parameter values are same as given in Eq. (5.2). Figure 8 corroborates that introduction of top predator even stabilizes an oscillatory system. Also, in both the above cases, introduction of top predator induces cascading effect. As top predator limits the density of mesopredator, the density of prey species increases. Phase portrait of systems (2.1) and (2.2) for k 1 = 1.9965 is plotted in Fig. 9. One interior equilibrium of system (2.1) is saddle and the other unstable one is surrounded by a stable (magenta color) and an unstable limit cycle (green color). Region inside the unstable limit cycle serves as the basin of attraction of the stable limit cycle, which is very narrow. Otherwise the trajectories finally accumulates at origin, and in the meantime the density of prey species falls below our considered x c . In such a situation, top predator introduction magnifies the basin of attraction of corresponding interior equilibrium many fold, mop up the oscillation, and drives the trajectories to stabilize at higher equilibrium density of prey species. Here also, some region of positive quadrant lies outside the basin of attraction of interior equilibrium, and trajectories starting in this region ultimately terminate at origin. We have already observed in Fig. 2(a), higher strength of fear triggers population collapse. So in Fig. 10, we consider k 1 = 2.2 and plot times series of system (2.1) up to t = 5 unit, when both the species rapidly move towards extinction. Introduction of top predator at that instance mediates the system to a stable state where all the species co-exist. All trajectories in Fig. 11(a) converge to trivial equilibrium for any initial start. In such a situation, species can be prevented from extinction, provided the initial point of the trajectories of system (2.1) lies inside the basin of attraction of stable equilibrium point (Fig. 11(b)). Thus, the density of a particular species in an ecosystem can be enriched by introducing a suitable top predator. Even the extinction of a species can be avoided by timely introducing top predator.
Thereafter, we plot the equilibrium curve for system (2.2) by varying fear parameter k 1 (Fig. 12). We observe that a stable branch exists up to large value of k 1 and then the system undergoes saddle-node bifurcation at k 1 = 14.08, where the stable branch collides with unstable branch and co-existence equilibrium vanishes. That means when perceiving predation risk, prey species deploy most of their effort 2450068-24 and time in anti-predation activities, growth rate becomes very low and prey population gradually dies out. This further implies that the whole system fall out and all the species become extinct. Figure 13(a) represents the equilibrium curve in α 2 -x plane, obtained by varying the predation rate of top predator (α 2 ), considering k 1 = 1 and rest of the parameter values are same as in Eq. (5.2). The system undergoes transcritical bifurcation at α 2 = 0, and the system remains stable for positive α . α 2 = 0 is just mimicking the situation that the dynamics is governed by system (2.1), and then the system achieves its stable equilibrium (x * , y * ) = (803.598, 321.439). Now, as the top predator is introduced, and its predation rate takes any positive value, density of prey species increases from its previous state. However, for k 1 = 4, the system shows different kinds of dynamics (Fig. 13(b)). System (2.2) remains stable for relatively higher values of α 2 . It undergoes Hopf bifurcation at α 2 = 0.029, which is subcritical in nature (l 1 = 4.734 × 10 -6 ), and saddle-node bifurcation at α 2 = 0.0257. Therefore, the system experiences population collapse or species extinction if α 2 decreases below 0.029. Since we consider k 1 = 4, we can see from Fig. 2(a) that system (2.1) already experiencing population collapse. Therefore, top predator is introduced at some instance of time to re-instate the species. However, if the top predators are found to be weak, failed to capture mesopredator at some considerable rate, then the target of enriching x-species population cannot be achieved. Furthermore, to get more intriguing effect of k 1 and α 2 on the dynamics of system (2.2), we plot bifurcation diagram in k 1 -α 2 parametric plane (Fig. 14). The red, blue and black color curve, respectively, representing the saddle-node, Hopf and homoclinic curve, all of which converge to BT-bifurcation point. Below the saddle-node curve, the system has no co-existence equilibrium. Between saddlenode and Hopf curve, although two co-existence equilibria exist, both are unstable. Between Hopf and homoclinic curve, one co-existence equilibrium gains stability
How can we avoid the extinction of any species naturally? through Hopf bifurcation with an unstable limit cycle surrounding the equilibrium.
In such case, if the initial point lies in the interior of unstable limit cycle, then trajectories converge to the stable co-existence equilibria, otherwise the trajectories fall out. Above the homoclinic curve, the unstable limit cycle vanishes with the occurrence of homoclinic bifurcation and the system has only a stable and an unstable co-existence equilibria. Biologically, the figure interprets that if the prey species is more fearful and show strong anti-predation behavior, then in order to protect ecosystem from collapsing, predation rate of introduced top predator should be high enough. Now, to understand the combining impact of both fear parameters on the density of targeted species (prey) population, we draw equilibrium curve of system (2.2) with respect to fear parameter k 2 , for three different values k 1 , represented by Fig. 15. The complete dynamics of system (2.2) with respect to k 1 and k 2 is captured in Fig. 16. Similar kind of dynamics as explained for Fig. 14 is observed in this case also. When the strength of mesopredator induced fear on prey population is weak, density of prey population maintained at higher level, whatever the strength of top predator induced fear might be. However, if the cost of fear on prey population is relatively high, then the strength of top predator induced fear should be comparatively high for the stable co-existence of the species in the system. Parameters should lie above the black curve (homoclinic curve) for stable co-existence of all the species with density of prey population at higher level. Therefore, if the targeted species is less responsive toward mesopredator induced fear, the purpose of elevating prey density would be successful by the introduction of top predator. However, for more fearful prey species, the motive can be achieved if the strength of top predator induced fear on mesopredator is comparatively higher.
Although increasing mesopredator induced fear factor (k 1 ) leads to shrinking of basin of attraction of the stable equilibrium of system (2.2), but increase in 2450068-27 5.2 5.4 5.6 5.8 6 6.2 6.4 6.6 6.8 0 0.005 0.01 0.015 0.02 Fig. 16. (Color online) Bifurcation diagram of system (2.2) in k 1k 2 parametric plane. Rest parameter values are same as given in Eq. (5.2). The red, blue and black color curves, respectively, representing the saddle-node, Hopf and homoclinic curve. 2450068-28 Int. J. Biomath. Downloaded from www.worldscientific.com by WEIZMANN INSTITUTE OF SCIENCE on 11/11/25. Re-use and distribution is strictly not permitted, except for Open Access articles. How can we avoid the extinction of any species naturally?
predation rate (α 2 ) and top predator induced fear (k 2 ) has a positive impact in expanding the basin of attraction.
Over the course of time, many species have gone extinct and a large number of species have been declared endangered. Sometimes, the species population in a habitat becomes low due to over-predation augmented with some environmental factors. Previously, the producers and the herbivores are believed to be main contributors in shaping an ecosystem and the ecologists conceived the concept of bottom-up approach that governs an ecosystem. With experimental evidences collected from various ecosystem, ecologists started to understand that an ecosystem is not totally controlled by the producers or primary consumers, those occupying the lower trophic level. The tertiary consumers or apex predators, those comprise the higher trophic levels of a food chain are also important drivers of ecosystem.
From the inception of Green W orld Hypothesis by Hairston et al. [20], the idea of top-down mechanism surfaced in the scientific community. From thereon, the role of predator in maintaining a healthy ecosystem is being studied more and more. In this study, we theoretically conceptualize the role of top predator in protecting a prey species from extinction, in presence of mesopredator. The main idea is that if the targeted species (prey) population in the considered habitat falls below a certain threshold, top predator would be introduced in that habitat. We can call it as species population enrichment in a habitat. For the proposed mathematical model, feasibility conditions of interior equilibrium and its stability conditions are obtained. Following are the results of our study obtained theoretically and numerically:
(1) For system (2.1), in the absence of top predator, the density of prey species is always decreasing with the increasing strength of fear and predation rate. Also, the basin of attraction of the stable co-existence equilibrium is observed to be shrinking critically with increasing fear parameter and predation rate. (2) Beyond certain threshold value of fear parameter or predation rate, population collapse and both species extinct. (3) Our theoretical study suggests that the elevation of density of lower trophic species in the considered region can be possible with the introduction of a top predator, which feed upon the mesopredator only but not on prey. (4) Also top predator introduction can eliminate the persistent periodic oscillation of species density and drive the system to stable state. (5) More importantly, species extinction can be prevented by timely introduction of top predator in the habitat. (6) If the impact of predation induced fear is very high on prey species, then in order to reinstate the density of prey species at higher level either the predation rate of top predator should be relatively high or the effect of fear on mesopredator is strong enough.