Chaotic dynamics of a tri-topic food chain model with Beddington–DeAngelis functional response in presence of fear effect

The most important fact in the field of theoretical ecology and evolutionary biology is the strategy of predation for predators and avoidance of prey from predator attack. A lot of experimental works suggest that the reduction of prey depends on both direct predation and fear of predation. We explore the impact of fear effect and mutual interference among predators into a three-species food chain model. In this manuscript, we have considered a tri-topic food chain model with Beddington–DeAngelis functional response between interacting species, incorporating the reduction of prey and intermediate predator growth due to the fear of intermediate and top predator, respectively. We have provided parametric conditions for existence of biologically feasible equilibria as well as their local and global stability. We have established conditions of transcritical, saddle-node and Hopf bifurcation about different equilibria. Finally, we have performed some numerical investigations to justify analytical findings.


Introduction
Analysis of dynamical activities of prey-predator interaction is one of the momentous themes for researchers in mathematical biology and theoretical ecology. Uniform existence and universal importance make the prey-predator interplay an attractive field of investigation. Prey-predator interplay becomes the focus of theoretical ecologist and experimental biologist from the last few decades [1,2]. Many mathematical models of prey-predator interactions have been formulated and investigated to identify the consumption and survival dynamics of species [3,4]. Prey-predator interactions may be governed by ordinary differential equations [3,5], fractional differential equations [6,7], partial differential equations [4,8], stochastic differential equations [9], delay differential equations [10] and difference equations also. Moreover, using mathematical models researchers are trying to identify different environmental effects such as Allee effect, refuge of species, harvesting, pattern formation, environmental fluctuations, etc.
In the early nineteenth century Malthus first formulated prey-predator interaction through mathematical modelling [11,12]. The celebrated Lotka-Volterra model was eventually enhanced by introducing logistic growth function for prey species [13,14], incorporating various functional response, environmental effects and these developments makes prey-predator interplay more and more realistic [15][16][17]. The joint impacts of inter-specific and intra-specific competition can affected the dynamics of competition outcomes [18]. The dynamics of prey-predator interaction can significantly be affected in presence of fear effect in the field of environmental biology and ecology [19]. The presence of predator populations has the most significant impact on a prey-predator system through direct predation as well as fear of predation. A large number of mathematical models have been concerned about considering the direct predation only [20][21][22][23][24] and reference therein. Those models are constructed with only prey dependent, both monotone and non-monotone functional responses [15,16,21,25,26], both prey-and predator-dependent functional response like Cowley-Martin type functional response [17], Beddington-DeAngelis type functional response [27][28][29], Monod-Haldane type functional response [30] and ratio-dependent functional response [21,31]. In 2020, Bai et al. explore the global stability of the boundary solution of a non-autonomous predator-prey system with Beddington-DeAngelis functional response [32]. In 2009, Takeuchi et al. studied the dissension along with investing time on taking care of juvenile and searching for food/nutrients of mature prey in absence of direct predation, while they assumed that matures accommodate their parental care time via learning [33]. Krivan (2007) investigated that exchange among foraging and predation based on classical Lotka-Volterra prey-predator system where either prey or predator or both species were adaptive to maximize their independent sturdiness [34].
The second one is that in presence of predator species, prey species may gently modify their behaviour because of the fear of predation threat. Appearance of predator may influence prey species more effectively than direct predation [35][36][37][38][39]. The birth rate of prey species has been subjected to modification through the introduction of fear effect [19,40]. A lot of field observation suggests to include the cost of fear in a prey-predator system, not only the direct predation [36][37][38][39]. It also may be the cause of anti-predator behaviour for prey species which plays a significant responsibility on demography of prey species [41]. For example, the reproduction physiology of elks (Cervus elaphus) influenced by wolves (Canis lupus) in Greate Yellowstone Ecosystem [42].
Prey species may transform their habitat from higher-risk zone to lower risk zone, to reduce the predation rate [38]. In 2011, Zanette et al. [39] observed during a whole breeding season that song sparrows (Melospiza melodia) reduces their reproduction due to fear of predators in their offspring season. This reduction is being occurred because of their anti-predator behaviour which persuades growth rate, in addition to their offspring endurance rates since female song sparrows laid a few eggs. Only some of those eggs can survive while most of the nestlings perished in a nest. They also observed that a variety of anti-predator behaviour is responsible for this effect. For example, frightened parents suckled their nestlings less, their nestlings were lighter and much more likely to perish. Correlational affirmation for birds [36,[43][44][45], Elk [42], Snowshoe hares [46] and dugongs [47] also display some indication that fear can interrupt prey-predator interactions. Very recently Elliott et al., [48] studied some field experiments on Drosophila melanogaster as prey and mantid as their predator species, to observe the effect of fear on populations robustness in relation to species density. They explored that in presence of mantid, the reproductive rate of drasophila reduces in both their breeding in addition to non-breeding seasons.
Depending on field experimental data, Wang et al., Zanette et al. [19,39] developed mathematical formulation for prey-predator interaction by introducing the cost of fear for prey because of predator species, where the level of fear shows a crucial role in the prey birth rate. They also observed that strong anti-predator behaviour or correspondingly most important cost of fear may reduce the risk of the existence of oscillatory behaviours and thus eliminate the scenarios "paradox of enrichment". They also displayed that cost of fear can stabilize the system by eliminating periodic behaviour as observed in prey-predator interactions. Also, periodic oscillations may arise, emerging from either subcritical or supercritical Hopf bifurcation under comparatively lower cost for fear [19]. Thus, the effect of fear can produce multi-stability in prey-predator interplays. In 2017, Wang and Zou [40] formulate a stage-structured prey-predator interaction with adap-tive avoidance for predator species by including fear of predator for prey species. They divided the prey population into juvenile and mature stage and constitute a system of delay differential system with maturation delay. Das and Samanta [49] analysed a stochastic preypredator system with the effect of fear due to predator on prey species and additional food is provided for the predator. Recently, Panday et al. [50] analysed a tritopic food chain model considering the effect of fear of top predator and intermediate predator on reproduction of intermediate predator and prey species, respectively. They studied the cost of fear on the stability dynamics of the model system.
In this manuscript, the Hastings-Powell model [51] has been modified incorporating the fear of intermediate and top predator on prey and intermediate predator, respectively, with Beddington-DeAngelis functional response for both species. We have established feasibility condition, local and global stability of different equilibrium point, one-dimensional bifurcation (transcritical, saddle-node, Hopf bifurcation) of the considered system and numerically explore the chaos control dynamics of the system. In some cases, the numerical results has been verified by using the Gillespie algorithm [52]. The organization of the manuscript is as follows: in the next segment, we have formulated the model system. The basic dynamical results such as positivity, boundedness property and persistence of model are provided in Sect. 3. Biologically feasible equilibria of the model system and parametric conditions of local and global stability are determined in Sect. 4. Sect. 5 is dedicated to studying transcritical bifurcation at an axial equilibrium point and conditions for the occurrence of saddle-node and Hopf bifurcation around coexistence equilibrium point along with stability direction of Hopf bifurcation. In Sect. 6, we perform some numerical simulations to justify our analytical findings, which also show the role of fear effect on the dynamics of prey-predator interactions. Finally, in Sect. 7, we summarize some biological indications from our analytical observation and possible future scope for upcoming research works.

Model formulation
In this research article, our main objective is to investigate the change of basic dynamical behaviour of a tri-topic food chain model in presence of fear due to the appearance of higher species. For this purpose, we start by considering the Hestings-Powell model [51]. In 1991 Hestings and Powell proposed a continuous tritopic food chain model by considering more naturalist logistic growth function for prey species and Holling type-II functional response in the form, X , Y and Z are population densities of prey, middle/intermediate predator and top predator population at any time T , respectively, and description of model parameters are given in Table 1.
In 2018, Pandey et al. [50], incorporated effect of fear in above considered tri-topic food chain model. They assumed that the intrinsic growth rate of prey population reduces due to the appearance of intermediate predator and the modified intrinsic growth rate of prey population is φ 1 is a monotonically decreasing function of both K 1 and Y . Similarly, the modified growth rate of intermediate is also a monotonically decreasing function of both K 2 and Z . Here, K 1 , K 2 are level of fear parameters of prey and intermediate predator population, respectively. On behalf of above consideration they investigate the dynamics of following modified tri-topic food chain model Fear functions φ 1 (K 1 , Y ) and φ 2 (K 2 , z) satisfies following ecological hypothesis : Level of fear parameter of intermediate predator Biomass Environmental protection for prey Dimensionless Intermediate predator mortality rate They discuss boundedness, positivity of system solutions and persistence of the considered model system.
Parametric conditions of local and global stability of different feasible equilibria are investigated by them also. However, in the study by Pandey et al. [50], it was assumed that prey-predator interplay depends only on prey density alone and it is of Holling type-II functional response. A huge number of detailed critical inspection on ecology and physiological evidence suggests that competition within predator species might be good for predator population under certain conditions in a deterministic environment [53]. To take account of the above observation, we have modified the model of Pandey et al. [50] by considering the Beddington-DeAngelis type function response, depending on both interacting species. Thus, the modified shape of above-considered model with Beddington-DeAngelis type function response is as follows: where B 1 , B 2 represents mutual interference among middle predators and top predators, respectively. To reduce number of system parameters, we introduce following dimensionless variables Then the above system reduces to following form: with the initial condition x(0) ≥ 0, y(0) ≥ 0, z(0) ≥ 0 and the dimensionless parameters are given in Table 1.
In next segment of the manuscript, we shall establish boundedness, positivity of system solutions and persistence of the system which will refer feasibility, well-posedness and existence of the system for a long time.

Positivity, boundedness and persistence of system solutions
In Theorem 1-4, we shall establish, respectively, positivity, boundedness of system solution and permanence of the model system (4).

Theorem 1
All system solutions of the model system (4) are positively invariant.
Proof The right hand side of the model system (4) are continuous functions of dependent variables x, y and z. After integration of equations of the system (4), we get It is apparent from above expressions that x(t), y(t) and z(t) remain non-negative for future infinite time if solution curve initiate from any interior point of R 3 Hence, R 3 + is positively invariant set for the system (4). Proof We consider the function L(t) as Differentiating the above relation with respect to t and using system (4), we get Therefore, for any arbitrary constant γ we have where γ = d ≤ mind 1 , d 2 ,(say) Using the theory of differential inequality, we have the following inequality Hence, as t → ∞, the above relation reduces to 0 < L ≤ M d . Therefore, all system solutions (x(t), y(t), z(t)) of the model system (4) with positive initial conditions i.e., initiate from R 3 + − {0} are uniformly bounded in the region: On the other hand, persistence of a system mathematically implies that minimum population densities of all three population are away from zero and bounded. Also, all populations can survive for a long time range. Dissipativeness: Biologically the dissipativeness refers upper boundedness of all population.

Theorem 3 The model system (4) is dissipative if
Proof From the system of Eq. (4), we have Applying standard comparison theorem [21], we can write from above inequality that lim t→∞ sup x(t) ≤ S 1 , lim t→∞ sup y(t) ≤ S 2 and lim t→∞ sup z(t) ≤ S 3 where S 1 , S 2 , S 3 are positive real solutions of the system of equations tively. Solving, the above system of equations we get . Here, S 1 , S 2 are positive real numbers and S 3 will be positive if conditions declared in the statement of the theorem is satisfied. Hence, the system is dissipative if Uniform Permanence : Biologically permanence means all three populations will survive in future time.
Mathematically it says that system solutions are always away from zero. (4) is permanent if c 1 > s 2 (1+k 1 s 2 ) and d 1 (a 1 s 1 +b 1 s 2 +c 1 ) < s 1 hold, where s 1 , s 2 are defined in the proof. Proof From the system of equations (4), we have

Theorem 4 The model system
Applying standard comparison theorem [21], we can write from above inequality that lim t→∞ inf x(t) ≥ s 1 , lim t→∞ inf y(t) ≥ s 2 and lim t→∞ inf z(t) ≥ s 3 , where (s 1 , s 2 , s 3 ) are positive real solution for (x, y, z) of the system of equations: Solving the above system of equations, we obtain s 3 is positive real root of the equation Here, s 2 is a positive real number and s 1 , s 3 will be positive real if c 1 > s 2 (1 + k 1 s 2 ) and d 1 (a 1 s 1 + b 1 s 2 + c 1 ) < s 1 . Hence, the system is permanent if conditions stated in the statement of the theorem are satisfied.

Equilibria and their stability analysis
Here, we find equilibria of the model system (4) with their existence criterion and studied their local and global stability also.

Equilibria of the model system
As we are studying an ecological system, we are mainly interested to feasible equilibrium points only. Feasible equilibria of the system (4) are summarized below: (1) The trivial equilibrium point E 0 (0, 0, 0) and it is always exists. (2) The predators-free/axial equilibrium point E 1 (1, 0, 0), which also always exists. x 2 satisfies the equation The above equation can have at most two positive real roots, hence the system (4) may have at most two planer equilibria.
where (x * , y * , z * ) is positive real solution of the system of equations Eliminating x, y from above equations we have verified z satisfies a polynomial equation of degree five. If z * be a positive real root of that polynomial equation then one can write It is obvious that if a 2 d 2 < 1 then y * will be positive for positive value of z * . In Fig. 1, we have presented the projection of various interior equilibria of the model system (4) in x −z plane with parametric restriction a 2 d 2 < 1 for different values of the parameter fear level of prey (k 1 ) with other parameter values are taken from Table 2. Here trivial, axial and one planer equilibrium points are always exist. The system have one interior equilibrium point for the value of k 1 = 0.02 (see Fig. 1a) and we represent it by E 31 . For k 1 = 0.09078075523 the system has two interior equilibria, namely E 31 and coincident interior equilibrium point E 3 (see Fig. 1b). Increasing the fear level to k 1 = 0.2, three interior equilibria arise (see Fig. 1c). The interior equilibrium point with higher, moderate and lower population density of top predator are represented by E 31 , E 32 and E 33 respectively. Then, the interior equilibrium point E 33 vanishes (see Fig. 1d for k 1 = 0.33) and remaining two interior equilibria exists. For k 1 = 0.3840366382 remaining two interior equilibria E 31 , E 32 coincide to the interior equilibrium point E 3 (see Fig. 1e). Figure 1f represents the existence of no interior equilibrium point for value of k 1 = 0.5. It is clear that the coincident interior equilibrium point of E 31 , E 32 and E 32 , E 33 are E 3 and E 3 respectively.

Local stability analysis of equilibria
The local stability of equilibria will be discussed by determining eigenvalues of the Jacobian matrix at corresponding equilibrium points. We set up stabilityinstability conditions of equilibria in the following theorems.  Table 2 Table 2 Empirical values of system parameters of the model (4)
Proof The Jacobian matrix of the model system (4) at the trivial equilibrium point E 0 is Then, eigenvalues of J (E 0 ) are 1, −d 1 and −d 2 . Since, one eigenvalue is positive and other two are negative, hence the trivial equilibrium point E 0 is always saddle. Solutions in the vicinity of this equilibrium point are unstable in nature.
Biologically this result is an evidence of the long survival of the system, i.e. all three species will not go to extinction concurrently.
Proof The Jacobian matrix of the model system (4) at thee axial equilibrium point E 1 is It is clear from expression of eigenvalues that the system will be locally asymptotically stable if Biologically the above result is most important because if the mortality rate of intermediate predator (say), then the intermediate predator will go to extinction and as a result the top predator will also go to extinction and only the prey population will survive in the system.

Theorem 7 The planer equilibria E
Proof The Jacobian matrix of model system (4) at the planer equilibria The characteristic equation of above Jacobian matrix is with ξ 1 = − (a 11 + a 22 ) and ξ 2 = a 11 a 22 − a 21 a 12 .
Then one eigenvalue of the Jacobian matrix J (E 2 ) is y 2 a 2 y 2 + c 2 −d 2 and this will be negative if y 2 < c 2 d 2 1−a 2 d 2 . Other two eigenvalues are negative or having negative real parts if ξ 1 > 0 and ξ 2 > 0.
Hence, the planer equilibria E 2 (x 2 , y 2 , 0) is locally asymptotically stable if the conditions stated in the theorem hold.
Proof The Jacobian matrix of model system (4) at inte- The characteristic equation of above Jacobian matrix is where According to the Routh-Hurwitz criterion [54], interior equilibria E * is locally asymptotically stable if conditions stated in the statement of the theorem are satisfied.
If C 1 C 2 = C 3 then the system losses its stability at interior equilibria E * (x * , y * , z * ) through Hopf bifurcation if the tranversility condition is satisfied.

Global stability analysis of equilibria
Here, we shall establish conditions and domain of global stability of different equilibria. The global stability of an equilibrium point implies that the point will follow the same nature irrespective of the initial population size ever, i.e. the solution curve will reach to the equilibrium point starting from anywhere in the state space.
Proof To establish the global stability of axial equilibrium point E 1 , we introduced the Lyapunov function V 1 as follows : Then, time derivative of the above expression is Therefore, 1 a 2 and the equality sign occurs when (x, y, z) = (1, 0, 0). Hence by using Lyapunov-Lasalle's invariance principle [55] we can say that the axial equilibrium point E 1 (1, 0, 0) is globally asymptotically stable.
Proof To establish the global stability of planer equilibria E 2 (x 2 , y 2 , 0), we construct the Lyapunov function V 2 as follows : Then, time derivative of the above expression along solution curve of the system (4) is 1 and at the planer equilibria Then, using these relations the above equation reduces to Hence by using Lyapunov-Lasalle's invariance principle [55] we can say that the planer equilibria E 2 (x 2 , y 2 , 0) is globally asymptotically stable.

Theorem 11 The interior equilibria E
Proof To established the global stability of interior equilibria E * (x * , y * , z * ) we introduced the Lyapunov function V * as follows : Then, time derivative of above expression is Using the Theorem 1 and Theorem 2, we consider that 2 4 . Hence, after some algebraic calculation equation (7) reduces to From the above expression we can easily observed that . Hence by using Lyapunov-Lasalle's invariance principle [55], we can conclude that the interior equilibria (x * , y * , z * ) is globally asymptotically stable if  Fig. 2 Solution trajectories of the model system (4) with various starting points considering The above condition of global stability of interior equilibria is too complicated in the sense of inferring any ecological justification. So, we will verify the global behaviour of the interior equilibria numerically. Considering values of system parameter as k 1 = 1, .01, we see that the system consist of non-interior equilibria E 0 (0, 0), E 1 (1, 0, 0), E 2 (0.11982, 0.21670, 0) and only one interior equilibrium point E * (0.78077, 0.13886, 1.21888). Among all equilibria, the only interior equilibrium point E * is stable and others are saddle or unstable spiral without any limit cycle. Then any solution trajectories with various starting point converges asymptotically to the stable interior equilibrium point E * (see Fig. 2). This implies the global asymptotic stability of the interior equilibrium point E * .

Local bifurcation analysis
In this segment, we shall now identify some changes in structural behaviour of the model system (4). The qualitative behaviour of system solution depending upon certain parameter changes as the parameter goes through a certain critical value. For changes in the value of parameter the vector field changes its structural behaviour and hence the system goes through a bifurcation. A certain value of the significant parameter at which the qualitative change of dynamics of our considered model system occurs is called the bifurcation point. Now, we discuss the occurrence of transcritical bifurcation at the axial equilibrium point, saddle-node bifurcation at the interior equilibrium point, Hopf bifurcation and the direction of Hopf bifurcation in vicinity of interior equilibria concerning some system parameter as the bifurcation parameter.

Transcritical bifurcation analysis
We have seen that axial equilibrium point E 1 (1, 0, 0 Proof We shall verify transversality conditions of the considered model system (4) at axial equilibrium point E 1 for the occurrence of transcritical bifurcation considering d 1 as bifurcation parameter. It is clear that the Jacobian matrix J (E 1 ) corresponding to axial equilibrium point E 1 has two negative eigenvalues and a zero eigenvalue at The Jacobian matrix J (E 1 ) and its transpose J (E 1 ) T has eigenvector V = (1, −(a 1 + c 1 ), 0) T and W = (0, 1, 0) T corresponding to zero eigenvalue, respectively. Now Hence, all conditions of Sotomayer's theorem are satisfied and the system (4) experiences transcritical bifurcation as the parameter d 1 crosses the critical value To represent transcritical bifurcation graphically, we consider k 1 = 0.5, c 1 = 0.42, d 1 = 0.75, d 2 = 0.62 and other parameter values as given in Table 2. Then, the system contains trivial, axial and one planer equilibrium point. Here, the planer equilibria is a stable spiral, while others are saddle. Then increasing the value of parameter d 1 to d 1 = 0.95, the planer equilibrium point loses its feasibility criterion and the axial equilibrium point becomes a stable node. Thus, the axial equilibrium point exchanges its stability with planer equilibrium point through the destruction of planer equilibrium point, i.e. the system experiences transcritical bifurcation as the parameter d 1 goes through a critical value d [T C] 1 = 0.92593 (see Fig. 3). Ecologically, the above result has high significance on system dynamics since there is a critical value of the intermediate predator normal death rate d 1 , above which it goes to extinction and as a result the top predator also goes to extinction. Thus, only the prey population survive in the system at its highest population density.

Saddle-node bifurcation analysis
Suppose that for k 1 = k [SN ] 1 two interior equilibria of the model system (4) coincide. Now we investigate parametric condition, under which there are two interior equilibria for and among them one is a saddle point and other is a node. Also there is no interior equilibrium point for 1 . Thus, two interior equilibria collided and then annihilated. At the coincident interior equilibrium pointĒ * (x * ,ȳ * ,z * ), the Jacobian matrix has a zero eigenvalue and consequently the equilibrium point becomes non-hyperbolic. Thus, the system can goes through saddle-node bifurcation at the coincident interior equilibrium point E * as the parameter k 1 passes through critical value k 1 = k [SN ] 1 if conditions of Sotomayor's theorem are satisfied [56].
Proof First we rewrite the consider model system (4) into the forṁ The Jacobian matrix of the system (4) at coincident interior equilibrium pointĒ * (x * ,ȳ * ,z * ) has characteristic equation of the form Here, expressions of η i , i = 1, 2, 3 are depend on k 1 . Hence, det (J (Ē * )) = −η 3 is a function of k 1 and we can control the sign of η i , i = 1, 2, 3 changing the values of parameter k 1 . We can obtain critical value for which η 2 < 0, η 3 = 0 and consequently eigenvalues of the Jacobian matrix J (Ē * ) will be zero, one positive, one negative in order to obtain saddlenode bifurcation at the coincident interior equilibrium pointĒ * (x * ,ȳ * ,z * ).
In Fig. 4 we have presented the one parametric bifurcation diagram with regards to k 1 for (a) a 1 = 0.  Table 2. Here, the green, red dot represents the stable, unstable equilibrium point, respectively    Table 2. It is clear from the figures that the system shows more complex dynamics in second case, i.e. when the mutual interference among the top predators is higher, environmental protection to the intermediate predator and its prey handling times are comparatively lower. In Fig. 4a one saddle-node bifurcation point and in Fig.  4b two saddle-node bifurcation point arises.
Considering values of system parameters as given in Table 2 with k 1 = 1.3, a 1 = 0.7, b 2 = 0.4, c 2 = 1.2 we have drawn one parametric bifurcation diagram of top predator with regard to k 1 (see Fig. 4a). Then, the system (4) contain two interior equilibrium points along with one trivial, one axial and one planer equilibrium point. Among the two interior equilibria, one is stable node and other is saddle. Then as we increase value of parameter k 1 to k [SN ] 1 = 1.4743353695 both interior equilibrium points coincide into one and finally diminishes from the system if we again increase the parameter k 1 (for detailed dynamics see Fig. 5). It is observed from Fig. 4b that the system (4) contains one interior equilibrium point for k 1 < 0.09078075523 and at k 1 = 0.09078075523 the system contains two coexistence equilibria among them one is coincident, other is the previous one. Then, the system (4) contains three interior equilibria (among the two new interior equilibria one is saddle and other is unstable spiral with stable limit cycle) for 0.09078075523 < k 1 < 0.3. The  Table 2. Here, the green, red dot represents the stable, unstable equilibrium point, respectively system contains two interior equilibria (one is saddle and other is unstable spiral with stable limit cycle) for 0.3 ≤ k 1 < 0.3840366382 and then both interior equilibria coincide at k 1 = 0.3840366382, finally diminishes from the system for k 1 > 0.3840366382. Here, non-interior equilibria are always exist. Thus, the system goes through saddle-node bifurcation two times with regards to the bifurcation parameter k 1 at critical values k 1 = 0.09078075523, k 1 = 0.3840366382. Phase diagrams of model system (4) before and after bifurcations values are similar as in Figs. 10(a, d) and 11(a, c).
It is clear from Fig. 4a that with the increase of fear parameter k 1 the co-existence possibility of all species decreases. In Fig. 5, we have drawn the phase portrait of the model system (4) for different values of the parameter k 1 in the left, right and at the critical value of parameter k 1 with other parameter values as considered above. From the ecological point of view, we can say that as the fear of intermediate predator to prey increases then the prey population decreases and as a result it creates food crisis for both the predators. Consequently, the growth of both predators decreases and finally top predator population goes to extinction and population density of prey and intermediate predator varies periodically with time.
It is clear from Fig. 4b that with the increase of b 2 the co-existence region of all species with respect to k 1 decreases. Since with the increase of mutual interference among top predators, it becomes more busy in their intra-specific competition and less fear of middle predator is sufficient for the decrease of population density in top predator.

Hopf bifurcation analysis
In Hopf bifurcation, a periodic oscillatory solution of considered model system appear or disappear throughout a local change in stability of an equilibrium point. Now, we discuss the occurrence of Hopf bifurcation at an interior equilibrium point E * of the model system (4) concerning some bifurcation parameter. We shall vary the parameter k 1 , so as to obtain a Hopf bifurcation about the interior equilibrium point E * (x * , y * , z * ).
The characteristic equation of the Jacobian matrix J (E * ) of model system (4) at the interior equilibrium point E * (x * , y * , z * ) is given in (6). Expressions of C i ; i = 1, 2, 3 and C 1 C 2 − C 3 are depends on k 1 if we fix other parameters at constants. Then, values as well as sign of C i ; i = 1, 2, 3 and C 1 C 2 − C 3 can be controlled by changing values of parameter k 1 . Thus, Hopf bifurcation of the system (4) at interior equilibrium point E * (x * , y * , z * ) can be expected for a certain value k 1 = k * 1 , when C i > 0; i = 1, 2, 3 and C 1 C 2 − C 3 = 0. Then, from characteristic equation (6), we have Then, solving the above equation, we have roots as λ 1,2 = ±i √ C 2 , and λ 3 = −C 1 . Thus for any k 1 ∈ (k * 1 − , k * 1 + ), where is a pre assign positive number, roots of above equation are in general of the form : λ 1,2 (k 1 ) = α 0 (k 1 ) + ±iβ 0 (k 1 ) and λ 3 (k 1 ) = −C 1 (k 1 ). To occur Hopf bifurcation the transversality condition d dk 1 (Reλ i (k 1 ))| k 1 =k * 1 = 0, i = 1, 2 must be satisfied. Since λ 1 is a root of equation (6), hence we must have Comparing real and imaginary parts, we have Taking derivative of above equations with respect to k 1 , we get where From equation (8) and (9), we get α 0 (k 1 and hence using results where C i denotes derivative of C i with respect to k 1 . Now, we summarise the above discussion in the form of following theorem : then the system (4) exhibit a Hopf bifurcation, leading to a family of periodic solutions that bifurcates from E * (x * , y * , z * ) in vicinity of k 1 = k * 1 , where C 1 , C 2 , C 3 are defined above.

Stability direction of Hopf bifurcating periodic solution
We find the direction of Hopf bifurcation by using the theorem of normal form. We also derive the stability property of bifurcating periodic solutions of the model system (4) as discuss in the literature [57]. Then, eigenvectors v 1 , v 3 of the Jacobian matrix J (E * ) corresponding to eigenvalues λ 1 = iω, λ 3 = −C 1 respectively at k 1 = k * 1 as given by Then, using the transformation x = x * + p 11 x 1 + p 12 y 1 + p 13 z 1 , y = y * + p 21 x 1 + p 22 y 1 + p 23 z 1 , z = z * + p 31 x 1 + p 32 y 1 + p 33 z 1 , the model system (4) reduces to where − (x * + p 11 x 1 + p 12 y 1 + p 13 z 1 )(y * + p 21 x 1 + p 22 y 1 + p 23 z 1 ) Here, it is clear that the system (10) has an equilibrium point at E(0, 0, 0). The Jacobian matrix of system (10) at the equilibrium point We calculate values of h 11 , h 02 , h 20 , H 101 , H 110 , H 21 , g 11 , g 20 , ω, w 20 , w 11 and h 21 by using following relations : We find out the constants w 11 and w 20 solving two equations : C 1 w 11 = −h 11 , (D − 2iω)w 20 = −h 20 . Also, Another discriminant quantities which will determine the direction of Hopf bifurcation are listed below : Here, the sign of σ determines the direction Hopfbifurcation. The direction of Hopf-bifurcation is supercritical (or sub-critical) if σ > 0 (or σ < 0); δ determines stability of periodic solutions according as if δ < 0 (or δ > 0) then solutions are stable (or unstable); and the period of bifurcating periodic solutions is determined by T according as the period increases (or decreases) if T > 0 (or T < 0). Considering values of system parameters as given in Table 2, the system contains three interior equilibria along with three unstable non-interior equilibria. Among three interior equilibria, one is a stable spiral and others are unstable spiral. With the increase of parameter k 1 , the stable interior equilibrium point looses its stability and becomes center for the critical value of parameter k 1 at k * 1 = 0.30198352. If we further increase the parameter k 1 , then the equilibrium point becomes an unstable spiral and any solution trajectory starting in the vicinity of the equilibrium point goes to a stable limit cycle around it. Thus, the system (4) goes through supercritical Hopf bifurcation with respect to the bifurcation parameter k 1 at critical value k 1 = k * 1 . In Fig. 6, we have drawn phase portrait of the model system (4) for k 1 < k * 1 and k 1 > k * 1 with other parameters as given in Table 2. We have evaluated the first Lyapunov number σ of Hopf bifurcation at Hopf bifurcation point and the value is σ = 0.7349. Thus from the above discussion, we can say that the system (4) goes through supercritical Hopf bifurcation concerning k 1 as the bifurcation parameter.
Again the system (4) goes through sub-critical Hopf bifurcation with respect to bifurcation parameter k 2 at critical value k * 2 = 1.073028. We have drawn phase portrait of model system (4) for k 2 < k * 2 and k 2 > k * 2 in Fig. 7 with other parameter values as given in Table  2.
Thus, we can say that both fear parameters has significant role on the stability criterion of the coexistence equilibrium point of the considered system. Hence, the coexistence of all species in the system with fixed population density can be regulated by the fear of intermediate and top predators.

Numerical illustration
In this segment of the paper, we shall carry out some numerical simulations to examine different complex dynamics (like stable, one, two periodic and higher periodic or chaotic solution, etc.) varying the model parameters of the considered model. As the fear of top predator and intermediate predator reduces the growth rate of intermediate predator and prey population, respectively, they will change the dynamics of the food chain model abruptly. For this purpose, we draw different bifurcation diagram, phase portrait and time series solution of all the three species of the model system (4). In some cases, the phase portrait or time series solution are verified by drawing the time series solution of all the three species of the system using the Gillespie algorithm.
(a) Dynamics of the system for different values of k 1 and k 2 : To study the effect of fear level of both predators on the dynamics of the proposed food chain model different bifurcation curves are drawn in the k 2 − k 1 parametric plane with values of other system parameter as given in Table 2 (see Fig. 8). To discuss the existence and stability dynamics of different interior equilibria we have drawn the saddlenode, Hopf bifurcation curve for interior equilibria and the creation or distraction of one interior equilibrium point in the k 2 − k 1 parametric plane considering other parameter values as given in Table 2 in Fig. 8. In Fig. 8, the blue line is the creation or distraction of one interior equilibrium point and red, magenta lines are the saddle-node, Hopf bifurcation curves for interior equilibria respectively. Thus, the k 2 −k 1 parametric plane is divided into seven distinct subregion (D 1 , D 2 , D 3 , D 4 , D 5 , D 6 and D 7 ) with different number or nature of interior equilibria in each subregion. The trivial, axial and one planer equilibrium points are always exist throughout the plane where the trivial, axial equilibrium points are saddle and the planer equilibrium point is an unstable spiral for the considered parameter values.
In the subregion D 1 only one interior equilibrium point exists and it is an unstable spiral. A stable limit cycle arises in the phase space and all solution curve starting from anywhere in the phase space goes to the stable limit cycle (see Fig. 9a, b). Then increasing k 2 as enter in subregion D 2 , two new interior  Table  2. Here, green dot represents the stable equilibrium point and red dots represent unstable equilibria  Table  2. Here, green dot represents the stable equilibrium point and red dots represent unstable equilibria  Table  2. Here the blue, red and magenta lines are, respectively, the creation or distraction of one interior equilibrium point, saddle-node and Hopf bifurcation curve for interior equilibria  Table 2. Here, the red dot represents the unstable equilibrium point equilibria (one is saddle and other is unstable spiral) arise. All the three interior equilibria are unstable spiral here and a stable limit cycle arises around one of them (see Fig. 9c, d). But if we increase k 1 instead of k 2 then in the subregion D 7 , no interior equilibrium point exist, i.e. the interior equilibrium point diminishes and only the non-interior equilibria exist for the system (see Fig. 11c, d).
Then if we enter into the subregion D 4 from D 2 in Fig. 8 by crossing the Hopf bifurcation curve, the number of interior equilibria remain same but the interior equilibrium point with highest top predator  Table 2. Here, the green, red dot represents the stable, unstable equilibrium point, respectively population density becomes a stable spiral. Here, the phase space divided into two distinct basin of attractor where solution curve starting from any basin either goes to the stable interior equilibrium point or goes to the stable limit cycle around the planer equilibrium point, i.e. the bi-stability behaviour of the model system (4) occurs here (see Fig. 10a, b, c).
From subregion D 4 , if we decrease the parameter k 1 and enter into the subregion D 3 , we observe that both unstable interior equilibria vanishes but the bi-stability nature remains unaltered as previ-  Table 2. Here, the red dot represents the unstable equilibrium point ous, i.e. the remaining stable interior equilibrium point and a stable limit cycle around the planer equilibrium point exist here (see Fig. 10d, e, f). But if increase the parameter k 1 and enter into the subregion D 5 , only one unstable interior equilibrium point (saddle one) vanishes and the remaining two interior equilibria has the same nature as in the subregion D 4 . Bi-stability also occurs here and the phase space divided into two distinct attractor (see Fig. 10g, h, i). The bi-stability nature of the system (4) has high significance from the biological point of view. It suggests that the survival of population depends not only on values of system parameter, it also depends on the initial population density. Next, if we decrease the parameter k 2 and enter into the subregion D 6 from D 5 by crossing the Hopf bifurcation curve, we observe that the stable interior equilibrium point become an unstable spiral with a stable limit cycle and the other one remain same. Thus, both interior equilibria are unstable and any solution curve goes to the stable limit cycle (see Fig. 11a, b).
From the ecological point of view, we can say that both parameters k 1 , k 2 play a key role on the existence of all system populations. Then, with the increase of intermediate predator fear to prey, the density of prey decreases and the co-existence equilibrium points losses its feasibility criterion. Consequently, the top predator population goes to extinction while prey and middle predator population density vary periodically with time. Again with the increase of fear parameter k 2 , a lesser number of intermediate predator are accessible in the system by a top predator. As a result, top predator abundance will decrease, the consumption of prey by the intermediate predator will decrease but prey will be available for the intermediate predator. In this situation again the density of intermediate predator will increase and consequently the top predator abundance will increase. Thus, the system will follow the oscillatory behaviour and survival possibility of all species in the system will increase.
Thus from the above discussion one can conclude that both fear parameters (consequently fear effects) can regulate the co-existence of all three species for considered values of other system parameter. If the fear of top predator to intermediate predator increases the system goes to stable situa-   Fig. 16 shows the chaotic behaviour of the system (4) for k 1 = 0.02, second row of Fig. 16 exhibits the two periodic limit cycle oscillation for k 1 = 1.0. Then first row of Fig. 17 shows one periodic limit cycle oscilla- Here, the green, red dot represents the stable, unstable equilibrium point, respectively tion for k 1 = 1.8 and second row of Fig. 17 shows stable behaviour for k 1 = 6.0. (IV) Chaos control dynamics of the control parameter k 2 : Similarly, the fear of top predator on intermediate predator k 2 can also control the chaotic dynamics of model system (4). With the increase of fear parameter k 2 , the system enters into stable situation from its chaotic state through periodhalving bifurcation (see Figs. 18 and 19). The system shows chaotic, higher periodic oscillations in k 2 < 0.0198; then two periodic oscillations in the interval (0.0198, 0.0593); limit cycle oscillations in the interval (0.0593, 0.3931), and stable behaviour for k 2 ≥ 0.3931. In Figs. 18 and 19, we have drawn the phase portrait and time series solution using the RK4 algorithm and Gillespie algorithm of all the three species of the system (4) for different values of the parameter k 2 . The first row of Fig. 18 shows the chaotic behaviour of the system (4) for k 2 = 0.001, second row of Fig. 18 exhibits the two periodic limit cycle oscillation for k 2 = 0.04. Then first row of Fig. 19 shows one periodic limit cycle oscillation for k 2 = 0.1 and second row of Fig. 19 shows stable behaviour for k 2 = 2.0.
Ecologically, fear parameters k 1 , k 2 are highly significant since these two parameters can control the chaotic dynamics of the system and turn into a stable one. With the increase of fear parameters k 1 , k 2 , growth rate of prey, intermediate preda-  Fig. 20b. It is obvious from the figure that both fear parameters has stabilizing effect on the system dynamics (see Fig. 20b).

Conclusion
In this investigation we have explored the dynamics of a tri-topic Hasting-Powell food chain model incorporating the fear of intermediate predator on growth function of prey, the fear of top predator on growth function of intermediate predator and replacing the prey dependent Holling type II functional response by the ratio-dependent Beddington-DeAngelis functional response. In previous sections, we have shown that the cost of fear can influence local and global dynamics of the system. The positivity, boundedness and persistence of the system solution implies that any feasible solution is always away from zero and remains bounded. We find all possible feasible equilibria and their feasibility condition in expressions of system parameters. Local asymptotical stability of different equilibria are determined. We also discuss the global stability of different non-trivial equilibria. The transcritical bifurcation at predator-free/axial equilibrium point concerning the middle predator death rate as bifurcation parameter and saddle-node, Hopf bifurcation at interior equilibrium point concerning the fear of intermediate predator on prey population as bifurcation parameter are investigated. We also find the stability direction of Hopf bifurcating periodic solution concerning the fear of intermediate predator on prey population as the bifurcation parameter. Finally, we justify theoretical findings and the significance of different Here, the green, red dot represents the stable, unstable equilibrium point, respectively system parameters numerically considering hypothetical values of system parameters.
Considering moderate values (see Table 2) of system parameters, we see that the system has three interior equilibria. Then with the increase of fear parameter k 1 , the system becomes unstable from the stable situation but with the increase of fear parameter k 2 the system enters into a stable situation from the unstable one. There exists a set of values of the parameters k 1 , k 2 for which bi-stability arises with one stable interior equilibrium point and one stable limit cycle around the top predator-free/planer equilibrium point. Thus, the co-existence of all species or extinction of top predator depends not only on parametric restriction, it also depends on initial population density.
Again if the handling time of top predator, environmental protection for intermediate predator are high and mutual interference among intermediate, top predators are low, then the system shows highly com-plex dynamics like chaotic or higher periodic state, two periodic, one periodic and stable spiral situation.
Panday et al. [50] modified the three species Hasting-Powell food web model considering fear effect in growth function of prey and intermediate predator.
They shows that the cost of fear can stabilize the system from a chaotic situation through the period halving bifurcation. They also shows that if the cost of fear increases, then top predator population goes to extinction. They claimed that the fear effect can regulate the stability of the system. The existence of prey and predator species depends on the cost of fear.
But, here we observe that not only fear parameters, mutual interference among intermediate and top predators can also regulate the system dynamics. Mutual interference among predators can also reduce the chaotic behaviour of the system through period halving bifurcation and the system becomes stable. Thus along with the increase of fear parameter and Here, the red dot represents the unstable equilibrium point mutual interference among predator the system goes to the stable situation from chaotic one. Hence we have more options to reduce the chaotic dynamics of the modified Hasting-Powell food chain model and we will feel more relax to consider values of system parameters for a stable solution of the system. It is important to note that most of the research work considering the fear effect are two species preypredator model and there are only a nominal number of work considering three species prey-predator model. Intermediate predator foraging behaviour suppressed by the top predator and thereby mitigates the impacts of overconsumption of intermediate predator on prey. We showed this scenario by considering a tri-topic food chain model with fear of higher tropic level species and mutual interference among predators. We observe that fear and mutual interference among predators can make the system stable from the chaotic situation. Thus, fear and mutual interference among predators enhances the persistence and stability of the tri-topic food web model. Therefore to conserve biodiversity and maintaining the ecosystems we can manipulate fear and mutual interference by artificial vocalization.  20 Stability-instability region in a b 1 − b 2 , b k 1 − k 2 parametric plane of the system (4). The system exhibits chaotic or higher periodic, two periodic, one periodic oscillation and sta-ble spiral behaviour in maroon, purple, red and green region, respectively. Both parameters has stabilizing effect on the system dynamics n 12 = 0, (a 2ȳ * + b 2z * + c 2 ) 3 n 31 = 0. p 1 = l 11 θ 2 1 + l 22 θ 2 2 + 2l 12 θ 1 θ 2 , p 2 = m 11 θ 2 1 + m 22 θ 2 2 + m 33 θ 2 3 + 2m 12 θ 1 θ 2 + 2m 23 θ 2 θ 3 + 2m 31 θ 3 θ 1 , p 3 = n 22 θ 2 2 + n 33 θ 2 3 + 2n 23 θ 2 θ 3 .