Controlling chaos in three species food chain model with fear effect

1.   Introduction

Chaotic behaviour exists in many natural systems, such as weather and climate. We can study this behaviour using a nonlinear mathematical model. Controlling chaos or stabilizing chaotic population dynamics becomes a challenge. Small manipulation in a chaotic system can control chaos, for example, stabilize an unstable periodic orbit. This manuscript mainly concerns with the chaotic dynamics of three species food chain model. Hastings and Powell ^[1] shows the appearance of chaos in a natural food chain or web models. Many researchers ^[2,3,4,5] have studied the Hastings and Powell model ^[1] by incorporating some biological factors such as alternative food, Allee effects, refugia, disease, fear.

In ecology, the predator-prey interaction happens at higher trophic levels and predators have an impact on prey populations which may be direct and indirect or both ^[6]. In the direct effect, the predator predates prey while in the indirect effect, the predator induces fear in the prey population and this fear can change the prey's behaviour ^[7]. The changing in the behaviour of prey includes foraging and reproduction ^[8,9]. Zanette et al. ^[9] have done fieldwork on song sparrows during the breeding season while eliminating the direct predation to observe the impact of fear. A reduction of 40% reproduction has been found in the number of offspring due to the fear of predation alone. Due to the fear effect, the prey shows a variety of antipredator responses, including different psychological changes, habitat changes, and foraging. For example, birds respond to their antipredator defences, by escaping from their homes ^[10] and mule deer reduces foraging time due to predation danger of mountain lions ^[11].

The proposed model of Wang et al. ^[12] shows a reduction in the production of prey due to fear of the predator. They observed that the high level of fear could stabilize the system (by excluding population oscillations). Further, they found that the fear effect could change population oscillations from a supercritical Hopf bifurcation to subcritical Hopf bifurcation, which induces multiple limit cycles. The impact of fear has been investigated on prey-predator systems with refuge and allee effect ^[13,14]. Fear effect also has been observed on discrete, delayed, diffusive, and eco-epidemiological models ^{[15,16,17,18,19,20]}. In the fear-induced prey-predator models, most of the prey-predator systems are studied with Holling type-Ⅱ functional response. Pal et al. ^[21] proposed a prey-predator model, where the interaction between prey and predator has been followed by Beddington-DeAngelis functional response. Based on their results, they conclude that the fear effect has both stabilizing and destabilizing effects.

Chaotic dynamics in tritrophic food chain model has recently been incorporated by Panday et al. ^[4]. For the first time, they assumed that the growth rate of prey and the middle predator is suppressed due to the cost of fear of middle and top predator, respectively. Further, they conclude that fear factors can control chaotic dynamics.

Our goal in this manuscript is to incorporate the fear factor in an earlier studied model ^[22,23] and since the system is chaotic so our effort would be to stabilize the dynamics of a proposed model system with the help of fear factors.

2.   Model formulation

Here we consider a chaotic model ^[22,23], describes the interaction of the species in the food chain of prey, specialist predator and top specialist predator. The nonlinear model system is described as follows:

Here, $u_1, u_2$ and $u_3$ are the respective population densities of prey, intermediate predator and top predator. $a_0, a_1, b_0, d_0$, $d_1, d_2, d_3, w_0, w_1, w_2, w_3$ and $c_3$ are positive constants. $a_0$ represents the intrinsic growth rate of prey $u_1$, $a_1$ is the mortality rate of the predator $u_2$ in the absence of $u_1$ only. The parameters $w_0, w_1, w_2$ are the maximum value which per capita growth rate can attain. $d_0, d_1, d_2,$ and $d_3$ are half saturation constants corresponding to Holling type-Ⅱ functional response, $b_0$ is the rate of competition among individuals of prey, parameter $c_3$ is the mortality rate of the top predator in the absence of intermediate predator, and $w_3$ is measure of its assimilation efficiency.

It is observed from the field experiments that the production reduces due to fear effects. The following assumptions are made to incorporate fear effect in model system (2.1):

(ⅰ) We assume that due to fear of top predator $u_3$, the growth rate of intermediate predator $u_2$ reduces. Therefore the modified growth rate change to $\frac{w_1}{1 + {f}_2 u_3}$ (a monotonically decreasing function of both ${f}_2$ and $u_3$).

(ⅱ) We also assume that the growth rate of prey reduces due to fear of intermediate predator and the growth rate change to $\frac{a_0}{1 + {f}_1 u_2}$, (a monotonically decreasing function of both ${f}_1$ and $u_2$).

Here $g_1({f}_1, u_2) = \frac{1}{1 + {f}_1 u_2}$ and $g_2({f}_2, u_3) = \frac{1}{1 + {f}_2 u_3}$ are fear functions, which account for the cost of anti-predator defence due to fear, where ${f}_1$ and ${f}_2$ represent the fear parameters (level of fear) of prey and intermediate predator, respectively. From the biological point of view, it is reasonable to assume that

and

By using above assumptions, the model system (2.1) takes the form:

We analyze the model system (2.4) with the following positive initial conditions

The rest of the paper is organised as follows: In section 3, the existence of equilibria is discussed. The local stability and bifurcation analysis are investigated in sections 4 and 5, respectively. Numerical simulation is discussed in section 6. Finally, the manuscript ends with a discussion in section 7.

3.   Existence of equilibria

In this section, the existence of non-negative equilibria are discussed and further, stability analysis of these equilibrium points are established. The following non-negative equilibrium points are obtained:

(ⅰ) Trivial equilibrium point $E_0 (0, 0, 0)$, corresponds to the total extinction of prey and predator species.

(ⅱ) Axial equilibrium point $E_1 ({a_0}/{b_0}, 0, 0)$, corresponds to the extinction of the predator species.

(ⅲ) Boundary equilibrium point $E_2 (\tilde{u_1}, \tilde{u_2}, 0)$, corresponds to the extinction of the top predator species, where

and $\tilde{u_2}$ is the positive root of the quadratic equation

The planar equilibrium $E_2 (\tilde{u_1}, \tilde{u_2}, 0)$ exists if the following conditions are satisfied

(ⅳ) We obtain positive equilibrium point $E^{*}(u_1^{*}, u_2^{*}, u_3^{*})$, by solving these system of nonlinear equations:

Solving the following Eq (3.3), we obtain,

value of $u_1^{*}$ is obtained by solving the equation

and $u_3^{*}$ is obtained from the equation

From Eqs (3.4), (3.5) and (3.6), it is clear that, for the existence of positive equilibrium point following conditions should be satisfied:

Now, in order to study the behavior of solution near the equilibrium points, we compute the Jacobian matrix of the model system (2.4) at any point $(u_1, u_2, u_3)$, which is given by

where,

4.   Local stability analysis

For local asymptotic stability, the solutions must approach to an equilibrium point under initial conditions close to the equilibrium point. Here, the criterion for local stability of the equilibria are obtained by linearizing the model system around the corresponding equilibrium point.

Theorem 1.

(a) The trivial equilibrium point $E_0 (0, 0, 0)$ is always unstable.

(b) Axial equilibrium point $E_1 (a_0/b_0, 0, 0)$ is locally asymptotically stable if $a_1 > \frac{a_0 w_1}{a_0 + b_0 d_1}$.

(c) The equilibrium point $E_2 (\tilde{u_1}, \tilde{u_2}, 0)$ is locally asymptotically stable if $b_0 > \frac{w_0 \tilde{u}_2}{(d_0 + \tilde{u}_1)^2}$ and $c_3 > \frac{w_3 \tilde{u_2}}{d_3 +\tilde{u_2}}$.

Proof.

(a) The Jacobian matrix at trivial equilibrium point $E_0 (0, 0, 0)$ is given by

and the eigenvalues of the Jacobian matrix $J_{E_0}$ are $a_0, -a_1$ and $-c_3$. So the equilibrium point ${E_0}(0, 0, 0)$ is always unstable.

(b) The Jacobian matrix at the axial equilibrium point $E_1 (a_0/b_0, 0, 0)$ is

and the eigenvalues of $J_{E_1}$ are $-a_0, -c_3$ and $-a_1 + \dfrac{a_0 w_1}{a_0 + b_0 d_1}$. Hence $E_1 (a_0/b_0, 0, 0)$ is local asymptotically stable if $a_1 > \dfrac{a_0 w_1}{a_0 + b_0 d_1}$.

(c) The Jacobian matrix associated with $E_2 (\tilde{u_1}, \tilde{u_2}, 0)$ is

where,

The corresponding characteristic equation of the matrix $J_{E_2}$ is

The characteristic Eq (4.1) have negative real parts if $P_5 < 0$ and $P_1 < 0$. Hence the equilibrium point $E_2 (\tilde{u_1}, \tilde{u_2}, 0)$ is locally asymptotically stable if the conditions stated in the Theorem 1(c) are satisfied.

Theorem 2. The positive equilibrium $E^{*}(u_1^{*}, u_2^{*}, u_3^{*})$ is locally asymptotically stable if the following conditions hold:

where, $\alpha, \beta$ and $\gamma$ are given in (4.6).

Proof. The Jacobian matrix associated with the positive equilibrium $E^{*}(u_1^{*}, u_2^{*}, u_3^{*})$ is

and corresponding characteristic equation is

where,

with

The interior equilibrium will be locally asymptotically stable if the coefficients of the characteristic Eq (4.5), $\rho_1, \rho_2$ and $\rho_3$, are satisfy the Routh-Hurwitz stability criterion, i.e. $\rho_1 > 0, \rho_3 > 0$ and $\rho_1\rho_2-\rho_3 > 0$.

Straightforward computation shows that $\rho_1 > 0$ and $\rho_3 > 0$, if

which is true from the condition (4.2).

Also since,

which is true, using the condition (4.3).

Where,

In this work, our primary attention is to stabilize the chaotic dynamics using fear parameters. We observe from the numerical simulation that the model system obtains stable dynamics after a limit cycle oscillation. Therefore the occurrence of Hopf bifurcation is obvious.

5.   Hopf-Bifurcation analysis

Hopf bifurcation is a point or a critical value of a parameter where the system stability changes and periodic solution arises. In the next theorem, we prove the existence of Hopf bifurcation, and for this, we choose fear parameter $f_1$ as a bifurcation parameter.

Theorem 3. When the fear parameter $f_1$ crosses a critical value $f^{*}_1$, then the system enters into Hopf-bifurcation around the positive equilibrium $E^{*}$ if the following conditions hold:

Proof. Since we have seen that the interior equilibrium is locally asymptotically stable and we know that the model system loses its stability due to change in some parameter value. Hence, we choose the fear parameter ${f_1}$ as the bifurcation parameter. If there exists a critical value ${f_1}^{*}$ such that

For ${f_1} = {f_1}^{*}$ the characteristic Eq (4.5) must be of the form

the above equation has three roots $-\rho_1({f_1}^{*}), i\sqrt{\rho_2({f_1}^{*})}$ and $-i \sqrt{\rho_2 ({f_1}^{*})}$. To show Hopf bifurcation occurs at ${f_1} = {f_1}^{*}$, we need to satisfy transversality condition

For all ${f_1}$, the roots are in general of the form

Now, we shall verify the transversality condition

Substituting $\sigma_j(f_1) = \mu ({f_1}) \pm i \nu ({f_1}),$ into (5.2) and calculating the derivative, we have

where

Notices that $\mu({f_1}^{*}) = 0, \nu({f_1}^{*}) = \sqrt{\rho_2({f_1}^{*})}$, then we have

Solving $\mu'({f_1}^{*})$ from Eq (5.3) we get,

If $[\rho_1 ({f_1}^{*}) \rho_2({f_1}^{*})]' \neq \rho_3 '({f_1}^{*})$ and

Thus the transversality conditions hold and hence Hopf-bifurcation occurs at ${f_1} = {f_1}^{*}$.

From numerical simulation, we observed that the Hopf bifurcation occurs with respect to the fear parameter $f_1$ at the critical value $f^{*}_1 = 0.02769$ around the equilibrium point (24.58, 10.77, 5.839) and also occurs with respect to the fear parameter $f_2$ at the critical value $f^{*}_2 = 0.5268$ around the equilibrium (35.24, 10.77, 0.8819). The parameters used in the calculation, for the bifurcation points, are fixed as given in Figure 1.

6.   Numerical simulation

The main goal of this manuscript is to control the chaotic dynamic of the model system (2.1). For this purpose, we have proposed a mathematical model (2.4) by incorporating a biological factor of fear. Model system (2.1) has been studied in literature by many authors ^[22,23]. The model system (2.1) exhibits chaotic attractor, which can be seen in Figure 1. When we increase the value of the intrinsic growth rate of prey $a_0$, the system shows a chaotic dynamic in the range $1.96 \leq a_0 \leq 2$. These changes can be seen from a bifurcation diagram plotted in Figure 2. From Figure 3, we observe that the model system enters into a chaotic dynamic from a stable focus via limit cycle oscillation and period-doubling. It is noticed that when ${f}_1 = {f}_2 = 0$ then our proposed model (2.4) is same as an earlier studied model ^[22,23,24].

Next, we studied the impact of fear on the dynamics of the model system (2.4) and observed that in the absence of fear (i.e., ${f}_1 = {f}_2 = 0$), the model system exhibits chaotic dynamics. Further, by varying ${f}_1$ and ${f}_2$, we investigate the effect of fear on the system dynamics. Firstly, we investigate the impact of fear of intermediate predator (${f}_1$) on prey while the fear of top predator on the intermediate predator is absent, i.e., ${f}_2 = 0$. We observe that the increase in fear effect ${f}_1$ makes the system dynamics (2.4) stable from chaotic, which can be seen by a bifurcation diagram plotted in Figure 4. The system remains chaotic for the low cost of fear in prey growth. However, an increased cost of fear (${f}_1 \geq 0.03$) stabilizes the system dynamics from chaotic to stable focus, see Figure 5. From Figure 6, it is observed that above the threshold value of ${f}_1 = 0.08$, the extinction of top predator population is possible, whereas prey and intermediate predator shows stable coexistence. Next, we consider the fear of top predator $({f}_2)$, which reduces the growth rate of an intermediate predator where the fear of intermediate predator is absent, i.e., ${f}_1 = 0.$ It is observed that an increase in the fear parameter ${f}_2$ makes system (2.4) stable from a chaotic dynamics, see Figures 7 and 8. The bifurcation diagram illustrates that the chaotic dynamics observed for the low cost of fear while the stable dynamics observed at a higher cost of fear (${f}_2 \geq 0.6$). We have also observed that the large value of ${f}_2 ({f}_2 = 500)$, top predator goes to be extinct, whereas prey and intermediate predator population remain stable (can be seen in Figure 9).

Next, our investigation is to observe the dynamics of the model system (2.4) in the presence of both fear factors (${f}_1$ and ${f}_2).$ From Figure 10, we observe that for small values of fear factors the system dynamics remain chaotic and the increased value of ${f}_1$ or ${f}_2$ or both changes system dynamics from chaos to stable focus. The chaotic or higher periodic oscillations change to period-doubling, period-doubling to limit cycle and limit cycle to stable focus, as we increase the cost of fear factors.

7.   Conclusions

In this work, we have analysed a three species food chain model (2.4) incorporating fear factors in the classical model (2.1). Next, analytically, we have proved the existence of biological feasible equilibrium points. Further, the local stability analysis of the proposed model system (2.4) is done corresponding to existing equilibria. Since the model system dynamics transits from the limit cycle to stable focus, hence there occurs Hopf bifurcation. Therefore we have also done Hopf bifurcation analysis, where we have proved that the system enters into Hopf bifurcation with the fear parameters $f_1$ and $f_2$ as bifurcation parameters.

The model system (2.4) exhibits chaotic dynamics, and a chaotic attractor is observed (Figure 1). From Figures 2 and 3, we observe that increase in the intrinsic growth rate of prey $(a_0)$ leads to chaotic dynamics which can be controlled by fear factors. Further, from Figures 4 and 5, it is observed that the increase in fear factor ${f}_1$ makes the system stable from chaotic dynamics (even in the absence of cost of fear of top predator in intermediate predator). Even the large cost of fear causes a population to become extinct, see Figure 6. Similarly, increase in fear factor ${f}_2$ makes system (2.4) stable, which has been shown in Figures 7 and 8. Also, the large cost of fear causes the top predator population to become extinct (Figure 9).

It should be noted that these chaotic dynamics hold for a small cost of fear. If we increase the cost of fear, then the system dynamics change from chaos to stable focus. The model system remains chaotic for the low cost of fear, whereas when we increase the value of ${f}_1$ or ${f}_2$ or both, then the system tends to stable dynamics (see, Figure 10). Therefore, we can conclude that the fear parameters control the chaotic dynamics in a food chain model.

Acknowledgments

The research work of the first author is supported by the Council of Scientific & Industrial Research (CSIR, India) under the file no: 09/1058(0006)/2016-EMR-I. The authors are very grateful to anonymous referees for reading their paper carefully and for several constructive remarks and suggestions.

Conflict of interest

The authors declare that there is no conflicts of interest in this paper.