Study of Fractional Order Tri-Tropic Prey-Predator Model with Fear Effect on Prey Population

Subrata Paul¹, Animesh Mahata^2*, Supriya Mukherjee³, Prakash Chandra Mali⁴, Banamali Roy⁵

¹Department of Mathematics, Arambagh Government Polytechnic, Arambagh, India.
²Mahadevnagar High School, Maheshtala, Kolkata, India.
³Department of Mathematics, Gurudas College, Narkeldanga, Kolkata, India.
⁴Department of Mathematics, Jadavpur University, Kolkata, India.
⁵Department of Mathematics, Bangabasi Evening College, Kolkata, India.

DOI: 10.4236/apm.2022.1211050                

Abstract

In this manuscript, we have studied a fractional-order tri-trophic model with the help of Caputo operator. The total population is divided into three parts, namely prey, intermediate predator and top predator. In addition, the predator fear impact on prey population is suggested in this paper. Existence and uniqueness along with non-negativity and boundedness of the model system have been investigated. We have studied the local stability at all equilibrium points. Also, we have discussed global stability and Hopf bifurcation of our suggested model at interior equilibrium point. The Adam-Bashforth-Moulton approach is utilized to approximate the solution to the proposed model. With the help of MATLAB, we were able to conduct graphical demonstrations and numerical simulations.

Keywords

Prey-Predator Model, Stability, Fear Effect, Hopf Bifurcation

1. Introduction

Due to its worldwide existence and domination, the progressive connection between predator and prey has been studied. Fear of predators on prey has been shown to have a significant impact on anti-predator defenses and lower prey reproduction in recent zoological studies on terrestrial vertebrates. Prey that is afraid forage less, which may lower fertility, and they live by famine [1] [2]. Furthermore, the fear effect has a negative impact on the physiological state of young prey, which has a negative impact on their adult survival. All organisms in taxonomy face the risk of predation, and they display a range of anti-predator behaviours, including as adjustments to habitat utilisation, feeding behaviour, alertness, and physiological changes [3] [4] [5]. However, there is clear evidence that predator fear is just as important as direct eating. Fear is induced by such behavioral changes, which can lead to stress-related physiology and have a poor impact on physical health [6] [7] [8], which can affect reproductive health and survival. As a result, predator fear has long-term implications for population dynamics and ecosystems. Fear of predators on prey causes people to spend more time being watchful instead of foraging and looking for better food and reduced habitats [9] [10] [11]. As a result, impact of fear on reproductive success and adult body mass growth are still present. Because people become more attentive when density declines, fear of predators may amplify Allee effects. Individual attentiveness rises when group size diminishes, according to field study on social animals. Predation is causing small groups of prey to become extinct [12]. In population demographics, group defense is a common phrase that describes an occurrence in which predation is reduced as prey’s capacity to protect or hide themselves increases when they are in large numbers or form a group.

The fractional order derivative, unlike the conventional derivative, has an essential trait known as the memory effect. For biological systems, fractional order derivatives are connected to the entire time domain, whereas integer order derivatives reflect a change or a specific attribute at a certain moment [13]. As a result, the fractional order derivative is more suited to simulating memory difficulties in biological systems [14]. Recent study areas including ecological modeling, epidemiology, financial mathematics, and the physical and mathematical sciences have seen an increase in popularity of the fractional order system. The difference between our conventional integer order system and fractional order systems are illustrated in a number of articles on fractional order platforms [15] - [20]. Some ideas have been offered to show our research into various stability criterions such as limit cycle, Hopf bifurcations, global stability, and so on. The effect of fear has been taken up recently in the study prey-predator models [21] [22]. The authors in [23] studied the global asymptotic stability and hopf bifurcation in a homogeneous diffusive predator-prey System with Holling type II functional response. Hopf bifurcation analysis of the repressilator model is proposed in [24]. Our major goal is to combine the fear effect with Caputo fractional order model.

1.1. Motivation and Novelties of the Article

In order to replicate real-world issues, several innovative fractional operators with various properties have been designed. Moreover, the integer derivative has a local identity, whereas the fractional derivative has a global character. Numerous varieties of fractional derivatives, both with and without singular kernels, are available today. Leibniz’s query from 1695 marks the beginning of the fractional derivative. The fractional derivative also improves in the improvement of the system’s consistency domain. We have the derivatives of Caputo, Riemann-Liouville, and Katugampola for singular kernels [25] [26]. There are two varieties of fractional derivatives without singular kernels: the Caputo-Fabrizio fractional derivative, which has an exponential kernel, and the Atangana-Baleanu fractional derivative, which has a Mittag-Leffler kernel [27]. Numerous academic articles, monographs, and novels have provided evidence to support this claim; for instance, [28] - [32]. Motivated by the abovementioned works and the advantages of Caputo derivatives, this study built a prey predator model with fear effect in Caputo sense. Because the Caputo derivative allows for the inclusion of conventional starting and boundary conditions in the derivation and because the derivative of a constant is zero, as opposed to the Riemann-Liouville fractional derivative, it is particularly helpful for describing real-world problems. Non-local operators, which may represent non-localities and certain memory effects, are typically better suited for such situations since they can account for power law, fading memory, and overlap effects.

The objective of the current work is:

● To study a tri-tropic prey predator model in Caputo environment and its stability analysis.

● Determination of all equilibrium points.

● Existence of Hopf bifurcation of the proposed model.

● Numerical Solution by Adam-Bashforth method.

1.2. Structure of the Paper

We discuss some important definitions and characteristics of fractional derivatives related to this article in Section 2. We have formulated a tri-tropic prey predator model with fear effect in Section 3. Section 4 consists of the existence, uniqueness, non-negativity, boundedness of solutions of the system. We have studied the local stability and global stability in Section 5. Influence of fear on population density in the model system is described in Section 6. In Section 7, we have analyzed Hopf bifurcation of the model system with respect to the fear effect parameter. In Section 8, we have used MATLAB (2018a) for numerical simulations in order to demonstrate the validity of our mathematical conclusions. Section 9 consists of the conclusion.

2. Preliminaries

The definitions and features of fractional derivatives that we offer the reader are both informative and practical.

Definition 1. [33] “The Caputo fractional derivative of order $0<\alpha \le 1$ for the function $u:C^n\left[0,\infty \right]\to ℝ$ is defined as

${}^{C}D{}_t^{\alpha }\left(u\left(t\right)\right)=\frac{1}{\Gamma \left(n-\alpha \right)}{\displaystyle {\int }_0^t\frac{1}{{\left(t-z\right)}^{\alpha +1-n}}\frac{{\text{d}}^n}{\text{d}{z}^n}u\left(z\right)\text{d}z}$,

where $C^n\left[0,\infty \right]$ is a n tines continuously differentiable function and the Gamma function is defined by Γ( ) such that $n-1<\alpha <n$ ”.

Theorem 1. [34] “If ${}^{C}D{}_t^{\alpha }u\left(t\right)$ is piecewise continuous, then $L\left({}^{C}D{}_t^{\alpha }\left(u\left(t\right)\right)\right)=z^{\alpha }L\left(u\left(t\right)\right)-{\displaystyle {\sum }_{i=0}^{l-1}{z}^{\alpha -i-1}{u}^{\left(i\right)}\left(0\right)}$, $l-1<\alpha \le l\in \mathbb{N}$ ,where the Laplace transform is denoted by $L\left(g\left(t\right)\right)$ ”.

Theorem 2. [35] “One-parametric and two-parametric Mittag-Leffler functions are described as follows: $E_{a_1}\left(z\right)={\displaystyle {\sum }_{i=0}^{\infty }\frac{z^i}{\Gamma \left(a_{1}i+1\right)}}$ and $E_{a_1,a_2}\left(z\right)={\displaystyle {\sum }_{i=0}^{\infty }\frac{z^i}{\Gamma \left(a_{1}i+a_2\right)}}$, where $a_1,a_2\in {ℝ}^{+}$ ”.

Lemma 1. [36] “Let $0<\alpha \le 1,u\left(t\right)\in C\left[p,q\right]$ and if ${}^{C}D{}_t^{\alpha }u\left(t\right)$ is continuous in $\left[p,q\right]$, then $u\left(x\right)=u\left(p\right)+\frac{1}{\Gamma \left(\alpha \right)}{\left(x-p\right)}^{\alpha }\cdot {}^{C}D{}_t^{\alpha }u\left(z\right)$, where $0\le z\le x$, $\forall x\in \left(p,q\right]$ ”.

Note 1. “If ${}^{C}D{}_t^{\alpha }u\left(t\right)\ge 0$ ( ${}^{C}D{}_t^{\alpha }u\left(t\right)\le 0$ ), $t\in \left(p,q\right)$, then $u\left(t\right)$ is a non-decreasing (non-increasing) function for $t\in \left[p,q\right]$ ”.

Lemma 2. “Let us consider the fractional order system as

${}^{C}D{}_t^{\alpha }\left(Y\left(t\right)\right)=\Psi \left(Y\right),Y_{t_0}=\left(y_{t_0}^1,y_{t_0}^2,\cdots ,y_{t_0}^n\right),y_{t_0}^j,j=1,2,\cdots ,n$,

with $0<\alpha <1$, $Y\left(t\right)=\left(y^1\left(t\right),y^2\left(t\right),\cdots ,y^n\left(t\right)\right)$ and $\Psi \left(Y\right):\left[t_0,\infty \right]\to {ℝ}^{n\times n}$. For calculate the equilibrium points, we have $\Psi \left(Y\right)=0$. These equilibrium points are locally asymptotically stable iff each eigenvalue ${\lambda }_j$ of the Jacobian matrix

$J\left(Y\right)=\frac{\partial \left({\text{Ψ}}_1,{\text{Ψ}}_2,\cdots ,{\text{Ψ}}_n\right)}{\partial \left(y^1,y^2,\cdots ,y^n\right)}$ calculated at the equilibrium points satisfies $\left|\arg \left({\lambda }_j\right)\right|>\frac{\alpha \pi }{2}$ ”.

Lemma 3. “Assume that $u\left(t\right)\in {ℝ}^{+}$ is a differentiable function. Then, for any $t>0$, ${}^{C}D{}_t^{\alpha }\left[u\left(t\right)-u^{*}-u^{*}\ln \frac{u\left(t\right)}{u^{*}}\right]\le \left(1-\frac{u^{*}}{u\left(t\right)}\right){}^{C}D{}_t^{\alpha }\left(u\left(t\right)\right)$, $u^{*}\in {ℝ}^{+}$, $\forall \alpha \in \left(0,1\right)$ ”.

3. Model Formulation

Step-1: Particular species never occupy the entire space in nature; they are either prey or predators of other species. The dynamical prey-predator model makes an assumption about the functional response, which is quantified by the amount of prey taken per predator per unit time. Let $u\left(t\right),v\left(t\right)$ and $w\left(t\right)$ represent the density of the prey population, intermediate predator, and top predator, respectively, at any timet. To formulate the model system, the following assumption is made:

1) In the absence of a predator, the prey population $u\left(t\right)$ increases logistically at an inherent growth rate of r, with k representing the environment’s carrying capacity.

2) According to Holling Type-I functional response, the intermediate predator $v\left(t\right)$ consumes the prey species of with the predation rate $a_1$. With the predation $b_1$, the top predator $w\left(t\right)$ also consumes the prey according to the law of mass action (Holling Type-I functional response).

3) The intermediate predator’s energy conversion coefficient and natural date rate, respectively are $a_2$ and $m_1$.

4) The energy conversion coefficient of predation and the natural date rate of the top predator, respectively are $b_2$ and $m_2$.

5) ${\gamma }_1$ is the intra-specific predation rate of top predator.

Kar and Ghosh [37] suggested the model system in the following way depending on these assumptions:

$\begin{array}{l}\frac{\text{d}u}{\text{d}t}=ru\left(1-\frac{u}{k}\right)-a_{1}uv-b_{1}uw,\\ \frac{\text{d}v}{\text{d}t}=a_{2}uv-m_{1}v,\\ \frac{\text{d}w}{\text{d}t}=b_{2}uw-m_{2}w-{\gamma }_1{w}^2.\end{array}$ (3.1)

Setp-2: We have now changed the model by include a fear element that the intermediate predator causes in the prey. Then, we write down the model system

with the fear factor $g\left(\rho ,v\right)=\frac{1}{1+\rho v}$ where $\rho$ represents the fear level as follows:

$\begin{array}{l}\frac{\text{d}u}{\text{d}t}=\frac{r}{1+\rho v}u\left(1-\frac{u}{k}\right)-a_{1}uv-b_{1}uw,\\ \frac{\text{d}v}{\text{d}t}=a_{2}uv-m_{1}v,\\ \frac{\text{d}w}{\text{d}t}=b_{2}uw-m_{2}w-{\gamma }_1{w}^2.\end{array}$ (3.2)

The function $g\left(\rho ,v\right)=\frac{1}{1+\rho v}$, due to fear (felt by prey) of predator, the birth

rate of prey species is reduced. In biological aspects of $\rho ,v,g\left(\rho ,v\right)$, it is appropriate to assume:

$\begin{array}{l}g\left(0,v\right)=1,g\left(\rho ,0\right)=1,\underset{\rho \to \infty }{\lim }g\left(\rho ,v\right)=0,\underset{v\to \infty }{\lim }g\left(\rho ,v\right)=0,\\ \frac{\partial g\left(\rho ,v\right)}{\partial \rho }<0,\frac{\partial g\left(\rho ,v\right)}{\partial v}<0.\end{array}$

Step-3: We have analyzed the suggested model utilizing the Caputo derivative of order $0<\alpha \le 1$.

$\begin{array}{l}{}^{C}D{}_t^{\alpha }u=\frac{\stackrel{\vee }{r}{}^{\alpha }}{1+\stackrel{\vee }{\rho }{}^{\alpha }v}u\left(1-\frac{u}{\stackrel{\vee }{k}{}^{\alpha }}\right)-\stackrel{\vee }{a_1}{}^{\alpha }uv-\stackrel{\vee }{b_1}{}^{\alpha }uw,\\ {}^{C}D{}_t^{\alpha }v=\stackrel{\vee }{a_2}{}^{\alpha }uv-\stackrel{\vee }{m_1}{}^{\alpha }v,\\ {}^{C}D{}_t^{\alpha }w=\stackrel{\vee }{b_2}{}^{\alpha }uw-\stackrel{\vee }{m_2}{}^{\alpha }w-\stackrel{\vee }{{\gamma }_1}{}^{\alpha }{w}^2.\end{array}$ (3.3)

Step-4: Now, we rewrite the parameters as follows for the sake of computation convenience: $\stackrel{\vee }{r}{}^{\alpha }=r$, $\stackrel{\vee }{\rho }{}^{\alpha }=\rho$, $\stackrel{\vee }{k}{}^{\alpha }=k$, $\stackrel{\vee }{a_1}{}^{\alpha }=a_1$, $\stackrel{\vee }{b_1}{}^{\alpha }=b_1$, $\stackrel{\vee }{a_2}{}^{\alpha }=a_2$, $\stackrel{\vee }{m_1}{}^{\alpha }=m_1$, $\stackrel{\vee }{b_2}{}^{\alpha }=b_2$, $\stackrel{\vee }{m_2}{}^{\alpha }=m_2$, $\stackrel{\vee }{{\gamma }_1}{}^{\alpha }={\gamma }_1$.

With beginning circumstances, the improved model system (3.3) may eventually be stated as follows:

$\begin{array}{l}{}^{C}D{}_t^{\alpha }u=\frac{r}{1+\rho v}u\left(1-\frac{u}{k}\right)-a_{1}uv-b_{1}uw,\\ {}^{C}D{}_t^{\alpha }v=a_{2}uv-m_{1}v,\\ {}^{C}D{}_t^{\alpha }w=b_{2}uw-m_{2}w-{\gamma }_1{w}^2.\end{array}$ (3.4)

4. Analysis of the Model

This section investigates the existence, uniqueness, non-negativity, and boundedness of the proposed model.

4.1. Existence and Uniqueness

Theorem 4.1.1. There exists a unique solution of the proposed model (3.4) for each non-negative initial condition.

Proof: We are seeking for a sufficient condition for the presence and uniqueness of the proposed model (3.4) solutions in the region $\Pi \times \left(0,T\right]$, where,

$\Pi =\left\{\left(u,v,w\right)\in {ℝ}^3:\max \Vert u\Vert ,\Vert v\Vert ,\Vert w\Vert \le M\right\}$.

The method employed in [38] is used. Consider a mapping $F\left(Y\right)=\left(F_1\left(Y\right),F_2\left(Y\right),F_3\left(Y\right)\right)$ where $Y=\left(u,v,w\right)$ and $\bar{Y}=\left(\bar{u},\bar{v},\bar{w}\right)$ :

$F_1\left(Y\right)=\frac{r}{1+\rho v}u\left(1-\frac{u}{k}\right)-a_{1}uv-b_{1}uw,$

$F_2\left(Y\right)=a_{2}uv-m_{1}v,$

$F_3\left(Y\right)=b_{2}uw-m_{2}w-{\gamma }_1{w}^2.$

For any $Y,\bar{Y}\in \Pi$ :

$\begin{array}{l}\Vert F\left(Y\right)-F\left(\bar{Y}\right)\Vert =\left|F_1\left(Y\right)-F_1\left(\bar{Y}\right)\right|+\left|F_2\left(Y\right)-F_2\left(\bar{Y}\right)\right|+\left|F_3\left(Y\right)-F_3\left(\bar{Y}\right)\right|\\ =\left|\frac{r}{1+\rho v}u\left(1-\frac{u}{k}\right)-a_{1}uv-b_{1}uw-\frac{r}{1+\rho \bar{v}}\bar{u}\left(1-\frac{\bar{u}}{k}\right)+a_1\bar{u}\bar{v}+b_1\bar{u}\bar{w}\right|\\ +\left|a_{2}uv-m_{1}v-a_2\bar{u}\bar{v}+m_1\bar{v}\right|+\left|b_{2}uw-m_{2}w-{\gamma }_1{w}^2-b_2\bar{u}\bar{w}+m_2\bar{w}+{\gamma }_1{\bar{w}}^2\right|\\ \le r\left|\frac{u}{1+\rho v}-\frac{\bar{u}}{1+\rho \bar{v}}\right|+\frac{r}{k}\left|u^2-{\bar{u}}^2\right|+\left(a_1+a_2\right)\left|uv-\bar{u}\bar{v}\right|\\ +\left(b_1+b_2\right)\left|uw-\bar{u}\bar{w}\right|+m_1\left|v-\bar{v}\right|+m_2\left|w-\bar{w}\right|+{\gamma }_1\left|w^2-{\bar{w}}^2\right|\end{array}$

$\begin{array}{l}=r\left|\frac{\left(u-\bar{u}\right)+\rho \left(u\bar{v}-v\bar{u}\right)}{\left(1+\rho v\right)\left(1+\rho \bar{v}\right)}\right|+\frac{r}{k}\left|u-\bar{u}\right|\left|u+u\right|+(a_1+a_2)\left|uv-\bar{u}\bar{v}\right|\\ +\left(b_1+b_2\right)\left|uw-\bar{u}\bar{w}\right|+m_1\left|v-\bar{v}\right|+m_2\left|w-\bar{w}\right|+{\gamma }_1\left|w-\bar{w}\right|\left|w+\bar{w}\right|\\ \le r\left|\frac{u-\bar{u}}{\left(1+\rho v\right)\left(1+\rho \bar{v}\right)}\right|+r\left|\frac{\rho \left(u\bar{v}-v\bar{u}\right)}{\left(1+\rho v\right)\left(1+\rho \bar{v}\right)}\right|+\frac{2Mr}{k}\left|u-\bar{u}\right|+\left(a_1+a_2\right)\left|uv-\bar{u}\bar{v}\right|\\ +\left(b_1+b_2\right)\left|uw-\bar{u}\bar{w}\right|+m_1\left|v-\bar{v}\right|+m_2\left|w-\bar{w}\right|+2M{\gamma }_1\left|w-\bar{w}\right|\end{array}$

$\begin{array}{l}\le r\left(1+\frac{2M}{k}\right)\left|u-\bar{u}\right|+\{r\left|\frac{\rho \left(u\bar{v}-v\bar{u}\right)}{\left(1+\rho v\right)\left(1+\rho \bar{v}\right)}\right|+\left(a_1+a_2\right)\left|uv-\bar{u}\bar{v}\right|\\ \begin{array}{c}{}_{}{}^{}\\ \end{array}+\left(b_1+b_2\right)\left|uw-\bar{u}\bar{w}\right|+m_1\}\left|v-\bar{v}\right|+\left(m_2+2M{\gamma }_1\right)\left|w-\bar{w}\right|\\ \le G_1\left|u-\bar{u}\right|+G_2\left|v-\bar{v}\right|+G_3\left|w-\bar{w}\right|\\ \le G\Vert Y-\bar{Y}\Vert .\end{array}$

where $G=\max \left\{G_1,G_2,G_3\right\}$ and $G_1=r\left(1+\frac{2M}{k}\right)$, $G_2=r\left|\frac{\rho \left(u\bar{v}-v\bar{u}\right)}{\left(1+\rho v\right)\left(1+\rho \bar{v}\right)}\right|+\left(a_1+a_2\right)\left|uv-\bar{u}\bar{v}\right|+\left(b_1+b_2\right)\left|uw-\bar{u}\bar{w}\right|+m_1$, $G_3=m_2+2M{\gamma }_1$.

As a result, $F\left(Y\right)$ fulfils the Lipschitz requirement. As a consequence, the fractional order system (3.4) exists and is unique.

4.2. Positivity and Boundedness of Proposed Model

Theorem 4.2.1. The model system’s solutions are all non-negative.

Proof: Let us assume $Z\left(t_0\right)=\left(u\left(t_0\right),v\left(t_0\right),w\left(t_0\right)\right)\in {\Omega }_{+}$ be the initial solution and ${\Omega }_{+}=\left\{\left(u,v,w\right)\in \Omega :u,v,w\in {ℝ}^{+}\right\}$, where ${ℝ}^{+}=ℝ\cup \left\{0\right\}$. Let us choose a constant $ℂ$, $t_0\le t<C$ such that

$\{\begin{array}{l}u\left(t\right)>0,\text{when}{t}_0\le t<C\\ u\left(ℂ\right)=0\\ u\left({ℂ}^{+}\right)<0\end{array}$

From system (3.4), we have ${}^{C}D{}_t^{\alpha }u\left(t\right)=0$ at $u\left(ℂ\right)=0$.

Using Lemma 1, we get $u\left({ℂ}^{+}\right)=0$, that contradicts $u\left({ℂ}^{+}\right)<0$. So $u\left(t\right)\ge 0$ $\forall t\in \left[t_0,\infty \right)$. In similar way we have $v\left(t\right)\ge 0,w\left(t\right)\ge 0$ $\forall t\in \left[t_0,\infty \right)$.

Theorem 4.2.2. The model system (3.4) has bounded solutions.

Proof: Let the function

$X\left(t\right)=u\left(t\right)+\frac{a_1}{a_2}v\left(t\right)+\frac{b_1}{b_2}w\left(t\right)$. (4.1)

Differentiating with respect to time on the above function, we have

$\begin{array}{l}{}^{C}D{}_t^{\alpha }X\left(t\right)={}^{C}D{}_t^{\alpha }u\left(t\right)+\frac{a_1}{a_2}{}^{C}D{}_t^{\alpha }v\left(t\right)+\frac{b_1}{b_2}{}^{C}D{}_t^{\alpha }w\left(t\right)\\ =\frac{r}{1+\rho v}u\left(1-\frac{u}{k}\right)-a_{1}uv-b_{1}uw+a_{1}uv-\frac{a_1{m}_1}{a_2}v+b_{1}uw-\frac{b_1{m}_2}{b_2}w-\frac{b_1{\gamma }_1}{a_2}{w}^2\\ =\frac{r}{1+\rho v}u\left(1-\frac{u}{k}\right)-\frac{a_1{m}_1}{a_2}v+\frac{a_1{b}_1}{a_2}uw-\frac{b_1{m}_2}{b_2}w-\frac{b_1{\gamma }_1}{b_2}{w}^2\\ \le -f\left(u+\frac{a_1}{a_2}v+\frac{b_1}{b_2}w\right)-\frac{r}{\left(1+\rho v\right)k}{\left(u-k\right)}^2-\frac{b_1{\gamma }_1}{a_2}{w}^2+\frac{r}{1+\rho v}k\end{array}$

where, $f=\min \left\{r,m_1,{\gamma }_1\right\}$.

Therefore, ${}^{C}D{}_t^{\alpha }X\left(t\right)+fX\left(t\right)\le R$, where, $R=\frac{r}{1+\rho v}k$.

With the help of Laplace transform, we have,

$p^{\alpha }L\left(X\left(t\right)\right)-p^{\alpha -1}X\left(0\right)+fL\left(X\left(t\right)\right)=\frac{R}{p}$

ð $L\left(X\left(t\right)\right)\left(p^{\alpha +1}+f\right)=p^{\alpha }X\left(0\right)+R$

ð $L\left(X\left(t\right)\right)=\frac{p^{\alpha }X\left(0\right)+R}{p^{\alpha +1}+f}=\frac{p^{\alpha }X\left(0\right)}{p^{\alpha +1}+f}+\frac{R}{p^{\alpha +1}+f}.$ (4.2)

Taking inverse Laplace transform, we have

$X\left(t\right)=X\left(0\right)E_{\alpha ,1}\left(-f{t}^{\alpha }\right)+R{t}^{\alpha }{E}_{\alpha ,\alpha +1}\left(-R{t}^{\alpha }\right).$

According to Mittag-Leffler function,

$E_{c,d}\left(z\right)=z{E}_{c,c+d}\left(z\right)+\frac{1}{\Gamma \left(d\right)}$.

Hence, $X\left(t\right)=\left(X\left(0\right)-\frac{R}{f}\right)E_{\alpha ,1}\left(-f{t}^{\alpha }\right)+\frac{R}{f}$.

Thus ${\lim }_{t\to \infty }SupX\left(t\right)\le \frac{R}{f}$. (4.3)

And hence the model (3.4) is bounded above by $\frac{R}{f}$.

As a result, all of the system’s (3.4) solutions will be bounded in ${ℝ}_{+}^3$.

5. Stability Analysis

In this part, we’ll identify all of system (3.4)’s trivial and non-trivial equilibrium points, as well as their existence conditions.

5.1. Equilibrium Points and Existence Criteria

There are five different types of equilibrium points in system (3.4).

1) $E_0\left(0,0,0\right)$ : Trivial equilibrium,

2) $E_1\left(k,0,0\right)$ : Axial equilibrium,

3) $E_2\left(u_1,v_1,0\right)$ : Planar equilibrium where,

$u_1=\frac{m_1}{a_2}$ and $v_1$ is to be obtained from the equation

$\frac{r}{1+\rho v_1}\left(1-\frac{u_1}{k}\right)-a_1{v}_1=0.$

$\Rightarrow \rho v_1^2+v_1-\frac{r}{a_1}\left(1-\frac{m_1}{a_2}\cdot \frac{1}{k}\right)=0.$

Hence $v_1$ is positive if $4\rho \frac{r}{a_1}\left(\frac{m_1}{a_2}\cdot \frac{1}{k}-1\right)<1$.

4) $E_3\left(u_2,0,w_2\right)$ : Planar equilibrium points where,

$u_2=\frac{m_1}{a_2}$, $w_2=\frac{m_2}{{\gamma }_1}\left[\frac{b_2}{m_2}\frac{m_1}{a_2}-1\right]$.

Hence $w_2$ is positive if $\frac{b_2}{m_2}\frac{m_1}{a_2}>1$.

5) The interior equilibrium is $E^{*}\left(u^{*},v^{*},w^{*}\right)$, here $\left(u^{*},v^{*},w^{*}\right)$ are the positive roots of the following system of equations

$\begin{array}{l}{E}^{*}\frac{r}{1+\rho v}\left(1-\frac{u}{k}\right)-a_{1}v-b_{1}w=0,\\ a_{2}u-m_1=0,\\ b_{2}u-m_2-{\gamma }_{1}w=0.\end{array}$ (5.1)

After solving the above equations, we get:

$u^{*}=\frac{m_1}{a_2},w^{*}=\frac{m_2}{{\gamma }_1}\left[\frac{b_2}{m_2}\frac{m_1}{a_2}-1\right].$ (5.2)

Hence $u^{*}$ and $w^{*}$ is positive if

$\frac{b_2}{m_2}\frac{m_1}{a_2}>1$. (5.3)

And $v^{*}$ is to be determined from the following equation

$\frac{r}{1+\rho v}\left(1-\frac{u}{k}\right)-a_{1}v-b_{1}w=0$. (5.4)

The Equation (5.4) has a positive root for

${\left(a_1+b_1\rho w^{*}\right)}^2>4\left[b_1{w}^{*}-r\left(1-\frac{u^{*}}{k}\right)\right]\left(a_1\rho \right).$

5.2. Local stability

The Jacobian matrix of the model (3.4) at $\left(u,v,w\right)$ can be represent as

$J\left(u,v,w\right)=\left[\begin{array}{ccc}A& B& C\\ D& E& O\\ F& O& G\end{array}\right]$.

where, $A=\frac{r}{1+\rho v}\left(1-\frac{2u}{k}\right)-a_{1}v-b_{1}w$, $B=-\frac{ru\left(1-\frac{u}{k}\right)\rho }{{\left(1+\rho v\right)}^2}-a_{1}u$,

$C=-b_{1}u,D=a_{2}v,E=a_{2}u-m_1,F=b_{2}w,G=b_{1}u-m_2-2{\gamma }_{1}w$

Theorem 5.2.1. The system (3.4) always exhibits unstable behavior at $E_0\left(0,0,0\right)$.

Proof: The eigenvalues of the Jacobian matrix are given by ${\lambda }_1=r$, ${\lambda }_2=-m_1$, ${\lambda }_3=-m_2$. It follows that $\left|\arg \left({\lambda }_1\right)\right|<\frac{\alpha \pi }{2}$ and $\left|\arg \left({\lambda }_j\right)\right|=\pi >\frac{\alpha \pi }{2}$ ( $j=1,2$ ). So, the system (3.4) is unstable at $E_0\left(0,0,0\right)$.

Theorem 5.2.2. The system (3.4) is locally stable at $E_1\left(k,0,0\right)$ if $k>\frac{m_1+m_2}{a_2+b_2}$.

Proof: The eigenvalue is ${\lambda }_1=-r$. To solve the equation

${\lambda }^2+\left(m_1+m_2-a_{2}k-b_{2}k\right)\lambda +\left[\left(m_1-a_{2}k\right)\left(m_2-b_{2}k\right)\right]=0$, (5.5)

we find another two eigenvalues.

Therefore, the condition of negative roots of Equation (5.5) is

$m_1+m_2<\left(a_2+b_2\right)k.$

So the system (3.4) is locally stable if $k>\frac{a_2+b_2}{a_2+b_2}$.

Theorem 5.2.3. The model system (3.4) is locally stable at $E_2\left(u_1,v_1,0\right)$ if $\left(b_1{m}_1\right)<\left(a_2{m}_2\right)$.

Proof: The matrix(J) at $E_2\left(u_1,v_1,0\right)$ can be written as:

$V\left(E_2\right)=\left[\begin{array}{ccc}\frac{r}{1+\rho v_1}\left(1-\frac{2u_1}{k}\right)-a_1{v}_1& -\frac{r{u}_1\left(1-\frac{u_1}{k}\right)\rho }{{\left(1+\rho v_1\right)}^2}-a_1{u}_1& -b_1\cdot \frac{m_1}{a_2}\\ a_2{v}_1& 0& 0\\ 0& 0& \frac{b_1\cdot m_1}{a_2}-m_2\end{array}\right]$.

One of the eigenvalue of Jacobian matrix is ${\lambda }_1=\left(\frac{b_1\cdot m_1}{a_2}-m_2\right)$ and the other two eigenvalues are obtain from of the equation ${\lambda }^2-\lambda \left(\frac{r}{1+\rho v_1}\left(1-\frac{2u_1}{k}\right)-a_1{v}_1\right)+\left(\frac{r{u}_1\left(1-\frac{u_1}{k}\right)\rho }{{\left(1+\rho v_1\right)}^2}+a_1{u}_1\right)\left(a_2{v}_1\right)=0$. So the system (3.4) is locally stable if $\left(b_1{m}_1\right)<\left(a_2{m}_2\right)$.

Theorem 5.2.4. The model system (3.4) is locally asymptotically stable at $E^{*}\left(u^{*},v^{*},w^{*}\right)$ if $A_i>0$, for $i=1,2,3$ and $A_1{A}_2-A_3>0$.

Proof: The characteristic equation of the Jacobian matrix at $E^{*}\left(u^{*},v^{*},w^{*}\right)$ is given by

${\lambda }^3+A_1{\lambda }^2+A_2\lambda +A_3=0$. (5.6)

where,

$A_1=\frac{r{u}^{*}}{\left(1+\rho v^{*}\right)k}+{\gamma }_1{w}^{*}$,

$A_2=\frac{r{u}^{*}}{\left(1+\rho v^{*}\right)k}\cdot {\gamma }_1{w}^{*}+b_1{u}^{*}{w}^{*}\cdot b_2+a_2{v}^{*}\left[\frac{r{u}^{*}\left(1-\frac{u^{*}}{k}\right)\rho }{{\left(1+\rho v^{*}\right)}^2}+a_1{u}^{*}\right]$,

$A_3=a_2{v}^{*}{w}^{*}{\gamma }_1\left[\frac{r{u}^{*}\left(1-\frac{u^{*}}{k}\right)\rho }{{\left(1+\rho v^{*}\right)}^2}+a_1{u}^{*}\right]$.

The roots of (5.6) are negative or have a negative real component, according to Rowth-Hurwitz criterion if $A_1>0,A_3>0$ and $A_1{A}_2-A_3>0$.

5.3. Global Stability

In this section we discuss the global stability of the system (3.4).

Theorem 5.3.1. The system (3.4) is globally asymptotically stable at $E^{*}\left(u^{*},v^{*},w^{*}\right)$ if

$\left(v-v^{*}\right)\left[\frac{r}{\left(1+\rho v\right)\left(1+\rho v^{*}\right)}-a_2\right]+\left[\frac{r\rho \left(u{v}^{*}-v{u}^{*}\right)}{k\left(1+\rho v\right)\left(1+\rho v^{*}\right)}-\left(w-w^{*}\right)b_2\right]\ge 0$.

Proof: Let the function L at $E^{*}\left(u^{*},v^{*},w^{*}\right)$ :

$\begin{array}{l}L=\left[\left(u-u^{*}\right)-u^{*}\log \left(\frac{u}{u^{*}}\right)\right]+\left[\left(v-v^{*}\right)-v^{*}\log \left(\frac{v}{v^{*}}\right)\right]\\ +\left[\left(w-w^{*}\right)-w^{*}\log \left(\frac{w}{w^{*}}\right)\right]\end{array}$ (5.7)

Taking time derivative of the above equation is

${}^{C}D{}_t^{\alpha }L=\frac{u-u^{*}}{u}{}^{C}D{}_t^{\alpha }u+\frac{v-v^{*}}{v}{}^{C}D{}_t^{\alpha }v+\frac{w-w^{*}}{w}{}^{C}D{}_t^{\alpha }w$. (5.8)

Using Lemma 3, then

$\begin{array}{l}{}^{C}D{}_t^{\alpha }L\le \frac{u-u^{*}}{u}\left[\frac{r}{1+\rho v}u\left(1-\frac{u}{k}\right)-a_{1}uv-b_{1}uw\right]+\frac{v-v^{*}}{v}\left[a_{2}uv-m_{1}v\right]\\ +\frac{w-w^{*}}{w}\left[b_{2}uw-m_{2}w-{\gamma }_1{w}^2\right]\\ \le \left(u-u^{*}\right)\left[\frac{r}{1+\rho v}\left(1-\frac{u}{k}\right)-a_{1}v-b_{1}w\right]+\left(v-v^{*}\right)\left[a_{2}u-m_1\right]\\ +\left(w-w^{*}\right)\left[b_{2}u-m_2-{\gamma }_{1}w\right]\\ =\left(u-u^{*}\right)b_{2}u-m_2-{\gamma }_{1}w\left[\frac{r}{1+\rho v}\left(1-\frac{u}{k}\right)-a_{1}v-b_{1}w\right]\end{array}$

$\begin{array}{l}+\left(v-v^{*}\right)\left[a_{2}u-m_1\right]+\left(w-w^{*}\right)\left[b_{2}u-m_2-{\gamma }_{1}w\right]\\ \le \left(u-u^{*}\right)\left[\left\{\frac{r}{1+\rho v}\left(1-\frac{u}{k}\right)-a_{1}v-b_{1}w\right\}-\left\{\frac{r}{1+\rho v^{*}}\left(1-\frac{u^{*}}{k}\right)-a_1{v}^{*}-b_1{w}^{*}\right\}\right]\\ +\left(v-v^{*}\right)\left[\left\{a_{2}u-m_1\right\}-\left\{a_2{u}^{*}-m_1\right\}\right]\\ +\left(w-w^{*}\right)\left[\left\{b_{2}u-m_2-{\gamma }_{1}w\right\}-\left\{b_2{u}^{*}-m_2-{\gamma }_1{w}^{*}\right\}\right]\\ -\frac{r{\left(u-u^{*}\right)}^2}{k\left(1+\rho v\right)\left(1+\rho v^{*}\right)}-\frac{r\left(u-u^{*}\right)\left(v-v^{*}\right)}{k\left(1+\rho v\right)\left(1+\rho v^{*}\right)}-\frac{r\left(u-u^{*}\right)\rho \left(u{v}^{*}-v{u}^{*}\right)}{k\left(1+\rho v\right)\left(1+\rho v^{*}\right)}\end{array}$

$\begin{array}{l}+a_2\left(u-u^{*}\right)\left(v-v^{*}\right)+b_2\left(u-u^{*}\right)\left(w-w^{*}\right)-{\gamma }_1{\left(w-w^{*}\right)}^2\\ \le -\frac{r{\left(u-u^{*}\right)}^2}{k\left(1+\rho v\right)\left(1+\rho v^{*}\right)}-{\gamma }_1{\left(w-w^{*}\right)}^2\end{array}$

$\begin{array}{l}-\left(u-u^{*}\right)\left(v-v^{*}\right)\left[\frac{r}{\left(1+\rho v\right)\left(1+\rho v^{*}\right)}-a_2\right]\\ -\left(u-u^{*}\right)\left[\frac{r\rho \left(u{v}^{*}-v{u}^{*}\right)}{k\left(1+\rho v\right)\left(1+\rho v^{*}\right)}-\left(w-w^{*}\right)b_2\right].\end{array}$

So ${}^{C}D{}_t^{\alpha }L\le 0$, if $\left(v-v^{*}\right)\left[\frac{r}{\left(1+\rho v\right)\left(1+\rho v^{*}\right)}-a_2\right]+\left[\frac{r\rho \left(u{v}^{*}-v{u}^{*}\right)}{k\left(1+\rho v\right)\left(1+\rho v^{*}\right)}-\left(w-w^{*}\right)b_2\right]\ge 0$.

6. Influence of Fear on Population Density in the Proposed Model

The main objective of this section is to investigate the impact of fear all populations. To look at this, we must first separate each component of the fear factor $\rho$ at $E^{*}\left(u^{*},v^{*},w^{*}\right)$.

Now $u^{*}=\frac{m_1}{a_2}$, $w^{*}=\frac{m_2}{{\gamma }_1}\left[\frac{b_2}{a_2}\frac{m_1}{m_2}-1\right]$. Hence $u^{*}$ and $w^{*}$ is positive if $\frac{b_2}{a_2}\frac{m_1}{m_2}>1$, and $v^{*}$ is the solution of the equation

$R{v}^{*}{}^2+S{v}^{*}+T=0$, (6.1)

where $R=\rho \left(a_1\right)$, $S=\left[a_1+\rho b_1\left\{\frac{m_2}{{\gamma }_1}\left(\frac{b_2}{a_2}\frac{m_1}{m_2}-1\right)\right\}\right]$ and $T=b_1\left\{\frac{m_2}{{\gamma }_1}\left(\frac{b_2}{a_2}\frac{m_1}{m_2}-1\right)\right\}-r\left(1-\frac{m_1}{a_2}\cdot \frac{1}{k}\right)$.

Differentiating (6.1) w.r. to $\rho$, we have $\frac{\text{d}{v}^{*}}{\text{d}\rho }=-\frac{v^{*}\left[L{v}^{*}+M\right]}{2R{v}^{*}+S}$, where

$L=a_1,M=b_1\left\{\frac{m_2}{{\gamma }_1}\left(\frac{b_2}{a_2}\frac{m_1}{m_2}-1\right)\right\}$

and $v^{*}=\frac{1}{2\rho }{a}_1\left(-a_1-\rho a_2{b}_1\left(b_2{m}_1-a_2{m}_2\right){\gamma }_1+\sqrt{X+Y}\right)$, with

$X=-4\rho a_1\left\{-\left(1-\frac{a_2{m}_1}{k}\right)r+a_2{b}_1\left(b_2{m}_1-a_2{m}_2\right){\gamma }_1\right\}$,

$Y={\left\{a_1+\rho a_2{b}_1\left(b_2{m}_1-a_2{m}_2\right){\gamma }_1\right\}}^2$.

As a result, the sign of $\frac{\text{d}{v}^{*}}{\text{d}\rho }$ is always negative, suggesting that as the fear level

$\rho$ rises, $v$ decreases. Differentiating equations (5.2) with respect to fear level $\rho$, we have

$\frac{\text{d}{u}^{*}}{\text{d}\rho }=0\text{and}\frac{\text{d}{w}^{*}}{\text{d}\rho }=0.$

As a result, we observed that as fear levels rise, the density of intermediate predators decreases.

7. Hopf-Bifurcation

Theorem 7.1. Let us consider ${}^{C}D{}_t^{\alpha }u=h\left(\rho ,u\right)$, where $0<\alpha <1,u\in {ℝ}^3$. A fractional order Hopf bifurcation is proposed in which the model (3.4) undergoes a Hopf bifurcation through $E^{*}$ at the value $\rho ={\rho }^{*}$ if

● The Jacobian matrix has two complex-conjugate eigenvalues ${\lambda }_{1,2}$.

● ${\theta }_{1,2}\left(\alpha ,{\rho }^{*}\right)=0$.

● ${\frac{\partial {\theta }_{1,2}}{\partial \rho }|}_{\rho ={\rho }^{*}}\ne 0$, where ${\theta }_i\left(\alpha ,{\rho }^{*}\right)=\frac{\alpha \pi }{2}-\left|\arg \left({\lambda }_i\left(\rho \right)\right)\right|,i=1,2$.

Proof: For $\rho ={\rho }^{*}$, the characteristic Equation (5.6) becomes

${\lambda }^3+{\phi }_1{\lambda }^2+{\phi }_2\lambda +{\phi }_1{\phi }_2=0$.

i.e., $\left({\lambda }^2+{\varphi }_2\right)\left(\lambda +{\phi }_1\right)=0$.

i.e., $\lambda =-{\varphi }_1,\pm i\sqrt{{\phi }_2}$.

Let ${\lambda }_1=i\sqrt{{\phi }_2},{\lambda }_2=-i\sqrt{{\phi }_2},{\lambda }_3=-{\phi }_1$.

For $\rho \in \left({\rho }^{*}-\delta ,{\rho }^{*}+\delta \right)$, where $\delta >0$, then we have

${\lambda }_1\left(\rho \right)={\sigma }_1\left(\rho \right)+i{\sigma }_2\left(k\right)$,

${\lambda }_2\left(\rho \right)={\sigma }_1\left(\rho \right)-i{\sigma }_2\left(k\right)$,

${\lambda }_3\left(\rho \right)=-{\phi }_1$.

Now the condition ${\frac{\partial {\theta }_{1,2}}{\partial \rho }|}_{\rho ={\rho }^{*}}\ne 0$ at $\rho ={\rho }^{*}$, is verify below.

Substitute ${\lambda }_1\left(\rho \right)={\sigma }_1\left(\rho \right)+i{\sigma }_2\left(\rho \right)$ into the characteristic equation and taking derivative with respect to $\rho$, we get

$A\left(\rho \right){{\sigma }^{\prime }}_1\left(\rho \right)-B\left(\rho \right){{\sigma }^{\prime }}_2\left(\rho \right)+S_1\left(\rho \right)=0$,

$B\left(\rho \right){{\sigma }^{\prime }}_1\left(\rho \right)+A\left(\rho \right){{\sigma }^{\prime }}_2\left(\rho \right)+S_2\left(\rho \right)=0$.

where

$A\left(\rho \right)=3{\sigma }_1^2\left(\rho \right)+2{\phi }_1\left(\rho \right){\sigma }_1\left(\rho \right)+{\phi }_2\left(\rho \right)-3{\sigma }_2^2\left(\rho \right)$,

$B\left(\rho \right)=6{\sigma }_1\left(\rho \right){\sigma }_2\left(\rho \right)+2{\phi }_1\left(\rho \right){\sigma }_2\left(\rho \right)$,

$S_1\left(\rho \right)={\sigma }_1^2\left(\rho \right){{\phi }^{\prime }}_1\left(\rho \right)+{{\phi }^{\prime }}_2\left(\rho \right){\sigma }_1\left(\rho \right)+{{\phi }^{\prime }}_3\left(\rho \right)-{{\phi }^{\prime }}_1\left(\rho \right){\sigma }_1^2\left(\rho \right)$,

$S_2\left(\rho \right)=2{\sigma }_1\left(\rho \right){\sigma }_2\left(\rho \right){{\phi }^{\prime }}_1\left(\rho \right)+{{\phi }^{\prime }}_2\left(\rho \right){\sigma }_2\left(\rho \right)$.

Noticing that, ${\sigma }_1\left({\rho }^{*}\right)=0,{\sigma }_2\left({\rho }^{*}\right)=\sqrt{{\phi }_2\left({\rho }^{*}\right)}$.

We have, $A\left({\rho }^{*}\right)=-2{\phi }_2\left({\rho }^{*}\right)$,

$B\left({\rho }^{*}\right)=2{\phi }_1\left({\rho }^{*}\right)\sqrt{{\phi }_2\left({\rho }^{*}\right)}$,

$S_1\left(\rho \right)={{\phi }^{\prime }}_3\left({\rho }^{*}\right)-{{\phi }^{\prime }}_1\left({\rho }^{*}\right){\phi }_2\left({\rho }^{*}\right)$,

$S_2\left(\rho \right)=-{{\phi }^{\prime }}_2\left({\rho }^{*}\right)\sqrt{{\phi }_2\left({\rho }^{*}\right)}$.

Therefore,

$\begin{array}{l}\frac{\text{d}}{\text{d}\rho }\left(Re\left(\lambda \left(\rho \right)\right)\right)\text{at}\rho ={\rho }^{*}\\ =-\frac{B\left({\rho }^{*}\right)S_2\left({\rho }^{*}\right)+A\left({\rho }^{*}\right)S_1\left({\rho }^{*}\right)}{{\left[A\left({\rho }^{*}\right)\right]}^2+{\left[B\left({\rho }^{*}\right)\right]}^2}\\ =-\frac{2{\phi }_1\left({\rho }^{*}\right){{\phi }^{\prime }}_2\left({\rho }^{*}\right){\phi }_2\left({\rho }^{*}\right)+2{{\phi }^{\prime }}_1\left({\rho }^{*}\right){\phi }_2\left({\rho }^{*}\right){\phi }_2\left({\rho }^{*}\right)-2{\phi }_2\left({\rho }^{*}\right){{\phi }^{\prime }}_3\left({\rho }^{*}\right)}{4\left[{\phi }_2^2\left({\rho }^{*}\right)+{\phi }_1^2\left({\rho }^{*}\right){\phi }_2\left({\rho }^{*}\right)\right]}\\ =-\frac{{\phi }_1\left({\rho }^{*}\right){{\phi }^{\prime }}_2\left({\rho }^{*}\right)-{{\phi }^{\prime }}_3\left({\rho }^{*}\right)+{{\phi }^{\prime }}_1\left({\rho }^{*}\right){\phi }_2\left({\rho }^{*}\right)}{2\left[{\phi }_2\left({\rho }^{*}\right)+{\left({\phi }_1\left({\rho }^{*}\right)\right)}^2\right]}\ne 0.\end{array}$

when, ${\phi }_1\left({\rho }^{*}\right){{\phi }^{\prime }}_2\left({\rho }^{*}\right)-{{\phi }^{\prime }}_3\left({\rho }^{*}\right)+{{\phi }^{\prime }}_1\left({\rho }^{*}\right){\phi }_2\left({\rho }^{*}\right)\ne 0$, and ${\lambda }_3\left({\rho }^{*}\right)=-{\phi }_1\left({\rho }^{*}\right)\ne 0$.

8. Numerical Discussions

To validate the theoretical findings presented in the earlier portions of this study, numerical simulations have been carried out using the modified Predictor-corrector approach [39] [40] in the Matlab framework. The Adams-Bashforth-Moulton scheme is employed for the numerical study.

${}^{C}D{}_t^{\alpha }{F}_j\left(t\right)=g_j\left(t,F_j\left(t\right)\right)$, $F_j^r\left(0\right)=F_{j0}^r$, (8.1)

$r=0,1,2,\cdots ,\lceil \alpha \rceil$, $j\in \mathbb{N}$

where $F_{j0}^r\in ℝ$, $\alpha >0$ and in the Caputo interpretation, $D_t^{\alpha }$ is equivalent to the popular Volterra integral equation.

$F_j\left(t\right)={\displaystyle {\sum }_{n=0}^{\lceil \alpha \rceil -1}{F}_{j0}^r\frac{t^n}{n!}}+\frac{1}{\Gamma \left(\alpha \right)}{\displaystyle {\int }_0^t{\left(t-u\right)}^{\alpha -1}{g}_j\left(u,F_j\left(u\right)\right)\text{d}u}$, $j\in \mathbb{N}$. (8.2)

We established model parameter values for numerical simulation purposes based on information from appropriate journal articles (see Table 1). This section is divided into four parts. The stability of our proposed model is discussed at $E_2$, $E_3$ and $E^{*}$ in Part 1. Part 2 delves into the dynamical behavior of all population of various fractional orders. Part 3 is to explore the Hopf bifurcation of the model system (3.4). Graph of mean density of all population with respect to the variation of $\alpha$ is discussed in Part 4.

Part 1:

The stability of our suggested model is discussed in this section. The parameter values used for the numerical simulations in Part 1 is provided in Table 1. Instead of the trivial $E_0$ and axial $E_1$, we are more concerned in the stability of the coexistence equilibrium $E_2$, $E_3$ and $E^{*}$. System (3.4) is shown in Figures 1-3 to be asymptotically stable for $\alpha =0.999$.

Part 2:

The parameter values in Table 1 are used to examine the dynamical behaviour

Table 1. Parameter values for numerical study.

(a)(b)

Figure 1. The system’s time series and phase portrait are consistent with Table 1 at $E_2$.

(a)(b)

Figure 2. The system’s time series and phase portrait are consistent with Table 1 at $E_3$.

(a)(b)

Figure 3. The system’s time series and phase portrait are consistent with Table 1 at $E^{*}$.

of the entire population. Figures 4(a)-(c) depict all populations’ behavior over time for different fractional orders α. Figure 4(a) depicts that the number of prey population increases when $\alpha$ changes from 0.85 to 0.95. We see in Figure 4(b) that number of intermediate predator population increases when α increases. Figure 4(c) depicts that the number of top predator population increases with time when α changes from 0.95 to 0.8.

Part-3:

In this part, it is explored whether the model system (3.4) exhibits a Hopf bifurcation with fractional order $\alpha =1$. The relevance of Hopf bifurcation is discussed using the following set of parametric variables.

The bifurcation analysis is investigated using the values of the parameters in Table 2. The model system’s unique endemic equilibrium $E^{*}\left(2.5,6.91,6.33\right)$ is obtained using the parameters in Table 2. Figure 5 shows the Hopf bifurcation diagram of the model system (3.4) with respect to the parameter $\rho$ taking $\alpha =1$. Also Figure 5 depicts that the model remains stable until $\rho$ crosses its threshold value $\rho ={\rho }^{*}=0.05$, when $\rho$ crosses its threshold value the model become unstable.

Part 4:

Figure 6 shows the change in mean density of prey, intermediate predator and top predator with respect to the variation of $\alpha$. Figure 6(a) depicts that mean density of prey population increases for $0\le \alpha \le 0.1$, decreases for $0.1<\alpha \le 0.8$ and oscillates for $\alpha >0.8$. Figure 6(b) shows that mean density of intermediate predator population increases for $0\le \alpha \le 0.45$, decreases for $0.45<\alpha \le 0.7$ and oscillates for $\alpha >0.7$. Figure 6(c) depicts that mean density of top predator population decreases for $0\le \alpha \le 0.8$.

(a)(b)(c)

Figure 4. Variation of $\alpha$ with respect to time of the model (3.4) corresponds to Table 1.

(a)(b)(c)

Figure 5. Bifurcation diagram of the system (3.4) corresponds with Table 2.

(a)(b)(c)

Figure 6. Graph of mean density of all population with respect to the variation of $\alpha$.

Table 2. Parameter values for study of Hopf bifurcation.

9. Conclusion

A fractional order prey-predator system with a fear impact on the prey population has been suggested. The non-negativity, boundedness, and uniqueness of the model system’s solutions have been demonstrated and investigated. The equilibrium points have been determined, and the stability of all of the model system’s possible equilibrium points has been studied both analytically and numerically. The system’s local and global stability are both established. Under constrained parametric circumstances, the system is found to be locally and globally asymptotically stable. The system is found to have Hopf bifurcations with respect to the parameter $\rho$. Local stability and Hopf bifurcation analysis have been the focus of our research. Fractional derivatives and integrals are a difficult idea to convey since they are derived from pure mathematics. The influence of predator anxiety on the population of prey is proposed in this research. A higher order indexing can be associated with weak memory, but a lower order classification can be associated with distant memory since the fractional order and memory are coupled. As a result, our research suggests that weak memories can help to improve the predator-prey system’s ability to cohabit peacefully, but powerful memories can actually make this situation worse. The discretization technique and FDE12 based on Adams-Bashforth-Moulton scheme are used to perform simulation studies. The concept of fractional calculus has nothing to do with any major geometrical meaning, such as function trend or convexity. Finally, we translate our mathematical findings into ecological terms as follows: a large number of prey refuges are produced in the system as a result of the prey species’ profound recollection of the exogenous effects that fear has on their life cycles. However, as we gradually lower the model system’s order, especially in the case of a low amount of predator-induced fear, the dynamics of the model system change away from its unstable behavior and toward stability. Consequently, our extensive mathematical findings show that weak memory might contribute to the stable existence of the predator-prey system whereas excessive memory of the species degrades the stability of the model system. In the suggested model (3.4), there is more work to be done, such as determining the species’ maturity to release harmful compounds into the environment. The model is now more realistic and intriguing as a result of these changes. We will leave this for further study.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Malthus, T.R. (1798) An Essay on the Principle of Population, and a Summary View of the Principle of Populations. Harmondsworth, Penguin.
[2]	Creel, S., Christianso, D., Lilley, S. and Winnie, J.A. (2007) Predation Risk Effects Reproductive Physiology and Demography of Elk. Science, 315, 960-960. https://doi.org/10.1126/science.1135918
[3]	Cresswell, W. (2011) Predation in Bird Populations. Journal of Ornithology, 152, 251-263. https://doi.org/10.1007/s10336-010-0638-1
[4]	Svennugsen, T.O., Holen, O.H. and Leimar, O. (2011) Inducible Defenses: Continuous Reaction Norms or Threshold Traits. The American Naturalist, 178, 397-410. https://doi.org/10.1086/661250
[5]	Preisser, E.L. and Bolnic, D.I. (2008) The Many Faces of Fear: Comparing the Pathways and Impacts on Non Consumptive Predator Effects on Prey Populations. PLOS ONE, 3, e2465. https://doi.org/10.1371/journal.pone.0002465
[6]	Sheriff, M.J., Krebs, C.J. and Boonstra, R. (2009) The Sensitive Hare: Sublethal Effects of Predator Stress on Reproduction in Snowshoe Hares. Journal of Animal Ecology, 78, 1249-1258. https://doi.org/10.1111/j.1365-2656.2009.01552.x
[7]	Elliott, K.H., Betini, G.S. and Norris, D.R. (2010) Experimental Evidence for Within- and Cross-Seasona Effects of Fear on Survival and Reproduction. Journal of Animal Ecology, 2016, 507-515. https://doi.org/10.1111/1365-2656.12487
[8]	Creel, S., Winnie, J.A. and Christianson, D. (2009) Gluco Corticoid Stress Hormones and the Effect of Predation Risk on Elk Reproduction. Proceedings of the National Academy of Sciences of the United States of America, 106, 12388-12393. https://doi.org/10.1073/pnas.0902235106
[9]	Creel, S. and Christianson, D. (2008) Relationships between Direct Predation and Risk Effects. Trends in Ecology & Evolution, 23, 194-201. https://doi.org/10.1016/j.tree.2007.12.004
[10]	Wirsing, A.J. and Ripple, W. (2010) A Comparison of Shark and Wolf Research Reveals Similar Behavioral Responses by Prey. Frontiers in Ecology and the Environment, 9, 335-341. https://doi.org/10.1890/090226
[11]	Schmitz, O.J., Beckerman, A.P. and O’Brien, K.M. (1997) Behavior Ally Mediated Trophic Cascades: Effects of Predation Risk on Food Web Interactions. Ecology, 78, 1388-1399. https://doi.org/10.1890/0012-9658(1997)078[1388:BMTCEO]2.0.CO;2
[12]	Manna, B., Paul, S., Mahata, A., Mukherjee, S. and Roy, B. (2022) Study of Prey-Predator Model Formulation and Stability Analysis. Advances in Intelligent Systems and Computing, 1422, 561-573. https://doi.org/10.1007/978-981-19-0182-9_57
[13]	Podlubny, I. (1999) Fractional Differential Equations. Academic Press, San Diego.
[14]	Mooring, M.S., Fitzpatrick, T.A., Nishihira, T.T. and Reisig, D.D. (2004) Vigilance, Predation Risk and the Allee Effect in Desert Bighorn Sheep. The Journal of Wildlife Management 68, 519-532. https://doi.org/10.2193/0022-541X(2004)068[0519:VPRATA]2.0.CO;2
[15]	Ahmed, E., El-Sayed, A. and El-Saka, H. (2007) Equilibrium Points, Stability and Numerical Solutions of Fractional-Order Predator-Prey and Rabies Models. Journal of Mathematical Analysis and Applications, 325, 5421-7553. https://doi.org/10.1016/j.jmaa.2006.01.087
[16]	Deshpande, A.S., Daftardar-Gejji, V. and Sukale, Y.V. (2017) On Hopf Bifurcational Dynamical Systems. Chaos, Solitons & Fractals, 98, 189-198. https://doi.org/10.1016/j.chaos.2017.03.034
[17]	Li, X. and Wu, R. (2014) Hopf Bifurcation Analysis of a New Commensurate Fractional-Order Hyper Chaotic System. Nonlinear Dynamics, 78, 279-288. https://doi.org/10.1007/s11071-014-1439-5
[18]	Das, M. Maity, A. and Samanta, G.P. (2018) Stability Analysis of a Prey Predator Fractional Order Model Incorporating Prey Refuge. Ecological Genetics and Genomics, 7, 33-46. https://doi.org/10.1016/j.egg.2018.05.001
[19]	Atangana, A. and Secer, A. (2013) A Note on Fractional Order Derivatives and Table of Fractional Derivatives of Some Special Functions. Abstract and Applied Analysis, 2013, Article ID: 279681. https://doi.org/10.1155/2013/279681
[20]	Li, H.L., Zhang, L., Hu, C., Jiang, Y.L. and Teng, Z. (2016) Dynamical Analysis of a Fractional-Order Predator-Prey Model Incorporating a Prey Refuge. Journal of Applied Mathematics and Computing 54, 435-449. https://doi.org/10.1007/s12190-016-1017-8
[21]	Wang, X., Zanette, L. and Zou, X.M. (2016) The Fear Effect in Predator-Prey Interactions. Journal of Mathematical Biology, 73, 1179-1204. https://doi.org/10.1007/s00285-016-0989-1
[22]	Zhang, H., Fu, S. and Wang, W. (2019) Impact of the Fear Effect in a Prey Predator Model Incorporating a Prey Refuge. Applied Mathematics and Computation, 356, 328-337. https://doi.org/10.1016/j.amc.2019.03.034
[23]	Wang, S., Yu, H., Dai, C. and Zhao, M. (2020) Global Asymptotic Stability and Hopf Bifurcation in a Homogeneous Diffusive Predator-Prey System with Holling Type II Functional Response. Applied Mathematics, 11, 389-406. https://doi.org/10.4236/am.2020.115028
[24]	Verdugo, A. (2018) Hopf Bifurcation Analysis of the Repressilator Model. American Journal of Computational Mathematics, 8, 137-152. https://doi.org/10.4236/ajcm.2018.82011
[25]	Caputo, M. and Fabrizio, M. (2015) A New Definition of Fractional Derivative without Singular Kernel. Progress in Fractional Differentiation and Applications, 1, 73-85.
[26]	Katugampola, U.N. (2011) New Approach to a Generalized Fractional Integral. Applied Mathematics and Computation, 218, 860-865. https://doi.org/10.1016/j.amc.2011.03.062
[27]	Atangana, A. and Baleanu, D. (2016) New Fractional Derivatives with Non-Local and Non-Singular Kernel Theory and Application to Heat Transfer Model. Thermal Science, 20, 763-769. https://doi.org/10.2298/TSCI160111018A
[28]	Wang, X., Zanette, L. and Zou, X. (2016) Modelling the Fear Effect in Predator-Prey Interactions. Journal of Mathematical Biology, 73, 1179-1204. https://doi.org/10.1007/s00285-016-0989-1
[29]	Pal, S., Pal, N., Samanta, S. and Chattopadhyay, J. (2019) Effect of Hunting Cooperation and Fear in a Predator-Prey Model. Ecological Complexity, 39, Article ID: 100770. https://doi.org/10.1016/j.ecocom.2019.100770
[30]	Sasmal, S.K. (2018) Population Dynamics with Multiple Allee Effects Induced by Fear Factors—A Mathematical Study on Prey-Predator Interactions. Applied Mathematical Modelling, 64, 1-14. https://doi.org/10.1016/j.apm.2018.07.021
[31]	Panday, P., Pal, N., Samanta, S. and Chattopadhyay, J. (2018) Stability and Bifurcation Analysis of a Three-Species Food Chain Model with Fear. International Journal of Bifurcation and Chaos, 28, 1-20. https://doi.org/10.1142/S0218127418500098
[32]	Wang, X. and Zou. X. (2017) Modeling the Fear Effect in Predator-Prey Interactions with Adaptive Avoidance of Predators. Bulletin of Mathematical Biology, 79, 1325-1359. https://doi.org/10.1007/s11538-017-0287-0
[33]	Kilbas, A., Srivastava, H. and Trujillo, J. (2006) Theory and Applications of Fractional Differential Equations. North-Holland Mathematics Studies, Vol. 204, Elsevier, Amsterdam.
[34]	Liang, S., Wu, R. and Chen, L. (2015) Laplace Transform of Fractional Order Differential Equations. Electronic Journal of Differential Equations, 2015, 1-15.
[35]	Kexue, L. and Jigen, P. (2011) Laplace Transform and Fractional Differential Equations. Applied Mathematics Letters, 24, 2019-2023. https://doi.org/10.1016/j.aml.2011.05.035
[36]	Petras, I. (2011) Fractional-Order Nonlinear Systems: Modeling Analysis and Simulation. Higher Education Press, Beijing. https://doi.org/10.1007/978-3-642-18101-6_3
[37]	Ghosh, B. and Kar, T.K. (2013) Possible Ecosystem Impacts of Applying Maximum Sustainable Yield Policy in Food Chain Models. Journal of Theoretical Biology, 329, 6-14. https://doi.org/10.1016/j.jtbi.2013.03.014
[38]	Li, H., et al. (2016) Dynamical Analysis of a Fractional-Order Predator-Prey Model Incorporating a Prey Refuge. Journal of Applied Mathematics and Computing, 54, 435-449. https://doi.org/10.1007/s12190-016-1017-8
[39]	Diethelm, K., Ford, N.J. and Freed, A.D. (2002) A Predictor-Corrector Approach for the Numerical Solution of Fractional Differential Equations. Nonlinear Dynamics, 29, 3-22.
[40]	Diethelm, K. (2003) Efficient Solution of Multi-Term Fractional Differential Equations Using P(EC)mE Methods. Computing, 71, 305-319. https://doi.org/10.1007/s00607-003-0033-3
[41]	Roy, J. and Alam, S. (2020) Fear Factor in a Prey-Predator System in Deterministic and Stochastic Environment. Physica A, 541, Article ID: 123359. https://doi.org/10.1016/j.physa.2019.123359
[42]	Barman, D., Roy, J. and Alam, S. (2020) Trade-Off between Fear Level Induced by Predator and Infection Rate among Prey Species. Journal of Applied Mathematics and Computing, 64, 635-663. https://doi.org/10.1007/s12190-020-01372-1