The hunting cooperation of a predator under two prey's competition and fear-effect in the prey-predator fractional-order model

1.   Introduction

In recent years, many studies have considered the ecological system's habitat from applied mathematics ^{[1,2,3,4,6,7,8,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,28,29,30,31,32,33,34,35,36,39,40,41,42]}. In these studeies, the fundamental research in mathematical modeling of the predator-prey model can be seen by Lotka and Volterra in ^[21,35], which opened many new aspects in modeling predator-prey interaction and showed the possibility of studying and predicting their dynamics.

In describing the habitat of the predator-prey model, functional response represents one of the essential terms. It shows the manner of interaction between prey and predator. More precisely, it expresses the attack methods and the quantity of the predator to the prey. In the literature, there are three forms of functional response: Type I, Type II, and Type III, where the predator consumption rate increases linear, hyperbolic or sigmoidal, respectively. Several mathematical forms represent this dynamics of the prey-predator model, which is $f(x) = ax$ of Holling Type I and $f(x) = ax/b+x$ of Holling Type II, where $a$, $b$ and $x$ are respectively the maximum predation rate, the half-saturation constant, and the prey biomass. On the other hand, the multi-species form of prey and their competition might increase the chance of the predator to attack and hunt them if the predator species itself is also in cooperation during hunting. Recently, Alves and Hilker ^[1] included a cooperation term to the rate of attack of the predator population and proposed the functional response $P(x, y) = (p+by)x$, where $x$ and $p > 0$ are respectively the prey density and the attack rate of the predator on the prey, while $b > 0$ describes the predator cooperation in hunting.

In this food chain cycle of the habitat, the fear of the prey started to become an important and realistic impact in the mathematical model to describe the predator-prey interaction, cooperation of the predators, and the fear of the prey. Many studies have been conducted to evaluate the effect of fear on prey population density ^[24,31,42] with different functions of response. However, it is new to study and examine the fear effect, the interaction of predator-prey, and the cooperation of predators in a single model.

Our following model considers an expanded form of Mukherjee ^[25] where we divide the ecological community into three compartments: $y(t)$ represents the predator, while there are two preys who are in competition and are denoted as $x_1(t)$ and $x_2(t)$, respectively. We consider a habitat where one predator compartment exists, which hunts $x_1(t)$ in cooperation, while species $x_2(t)$ is hunted individually. Also, we assume that both preys $x_1(t)$ and $x_2(t)$ are in competition which allows the predator to attack and hunt them easier. The collective hunting of predator $y(t)$ to species $x_1(t)$ activates a fear effect of the prey towards the hunter.

Therefore, this scenario is formulated as an ODE system as follows:

We equipped the system (1.1) with the initial condition

All the biological description of the parameters are given as below:

The modeling of prey-predator systems through fractional-order differential equations has many advantages. The nonlocal property of fractional-order models not only depends on the current state but also depends on its prior historical states. The transformation of an integer-order model into a fractional-order model needs to be precise with respect to the order of differentiation $\alpha$. However, a small change in $\alpha$ may cause a big change in the behavior of the solutions ^[41]. Fractional-order differential equations can model complex biological phenomena with non-linear behavior and long-term memory, which cannot be represented mathematically by integer-order differential equations (IDEs). For example, Bozkurt established the glioblastoma multiform (GBM)-immune system (IS) interaction using a fractional order differential equation system to include the delay time (memory effect) in ^[5]. Besides this, in many papers, researchers have proven the importance of fractional derivatives in modeling biological phenomena such as freedom towards ordering the derivative, dealing with species memory that has been achieved during their cycle life, genetic characteristics, and others (see for example ^{[4,6,8,9,13,27,30,37,38,41]}). Based on this information, we modify model (1.1) in a fractional-order form to present the study in a more natural state that suits the food chain life cycle requirements. Thus, we will consider the fractional-order ecological model such as

2.   Preliminaries

In this part, we want to introduce some fundamental properties of fractional-order differential equations, which will be helpful in the main sections.

Definition 2.1. ^[7] Caputo fractional derivatives can be given as follows

where $D^\alpha$ is Caputo differential operator of order $\alpha$, $d$ is the least integer, which is not less than $\alpha$, and $I^{\theta}$ is the Riemann-Liouville integral operator of order $\theta$, which is given by

where $\Gamma(\theta)$ is the Euler's Gamma function.

Theorem 2.2. ^[7] Consider the $N-$dimensional system

where $B$ is arbitrary constant $N\times N$ matrix and $\alpha\in(0, 1)$.

1) $x = 0$ is asymptotically stable if and only if all eigenvalues $\lambda_i, \; i = 1, 2, \ldots, N$ of $B$ satisfy $|arg(\lambda_i)| > \frac{\alpha\pi}{2}$.

2) $x = 0$ is stable if and only if all the eigenvalues of $B$ satisfy $|arg(\lambda_i)|\geq\frac{\alpha\pi}{2}$ and eigenvalues with $|arg(\lambda_i)| = \frac{\alpha\pi}{2}$ have same geometric multiplicity and algebraic multiplicity.

Theorem 2.3. ^[7] Consider the fractional order system

An equilibrium point $x^{*}$ is locally asymptotically stable if all the eigenvalues $\lambda_i$ of $J = \frac{\partial f(x^{*})}{\partial x}$ satisfy $|arg(\lambda_i)| > \frac{\alpha\pi}{2}$.

3.   Non-negativity and boundedness of solution in system (1.2)

This part investigates the existence and uniqueness of the solution in (1.2) and the non-negativity and boundedness.

Theorem 3.1. System $(1.2)$ shows a unique solution for non-negative initial conditions.

Proof. Let $E = \left\lbrace (x_1, x_2, y)\in R^3;max\{|x_1 |, |x_2 |, |y|\le K\}\right\rbrace$. In this proof, we use Hong et al. approach's ^[13]. Let $X = (x_1, x_2, y)$ and $H(X) = (H_1 (X), H_2 (X), H_3 (X))$, where we have

For $X, \bar{X} \in E$, it obvious to see from (3.1) that we get

where

From the analysis mentioned above, the Lipschitz condition of $M$ satisfies. Thus, we can conclude that the solution of system (1.2) exists and shows a unique solution.

Theorem 3.2. The boundedness and non-negative behavior of all solutions of system $(1.2)$ are permanent in $R_+^3$.

Proof. Let $L$ be the function defined as $L(t) = x_1 (t)+x_2 (t)+\frac{1}{c} y(t)$. Thus, for $\tau = min\{d \eta r_1-\kappa(1-\eta), m\}$, and $sup_{t\to \infty} x_2 (t) = r_2/a_2$, we have

By using the comparison theorem in ^[6], we obtain

where $G_\alpha$ is the Mittag-Leffler function. Using Lemma 5 and Corollary 6 of ^[6], we have

Hence, it is proven that all solutions of model (1.2) initiating in $R_+^3$ are in the region $\Gamma$, where

Now, we can see also the non-negativity of all solutions in system (1.2). From the first equation of (1.2) and $\Gamma$ we obtain,

where

By using the standard comparison theorem for fractional-order in Choi et al. ^[6] and the fact that $G_{\alpha, 1} (t) > 0$ for any $\alpha \in (0, 1)$ in ^[37], it follows that

Similarly, we can get from the second and third equations of (1.2) the following

Hence, obtain also $y(t) \ge 0$ and $x_2 (t) \ge 0$.

To analyse the stability conditions for the next section, we need to introduce the Jocabian matrix $J(x_1, x_2, y)$ associated to system (1.2), which is given by the following:

4.   Equilibria and their stability

System (1.2) has four equilibria:

1) Trivial equilibrium $P_0 (0, 0, 0)$, which exists always.

2) Prey equilibrium $P_1 (\frac{d-1}{a}, 0, 0)$, which exists provided that $d > 1$.

3) Predator-free equilibrium $P_2 (\bar{x}_1, \bar{x}_2, 0)$ where $\bar{x}_1 = \frac{1}{\beta_2} (r_2 -a_2 \bar{x}_2)$ and $\bar{x}_2 = \frac{a_1r_1+\beta_2d-\beta_2 r_1}{a_1a_2-\beta_1 \beta_2}$, which exists provided that: $a_1r_1+\beta_2d > \beta_2r_1$, $a_1a_2 > \beta_1 \beta_2$ and $r_2 > a_2 \bar{x}_2$.

4) The positive equilibrium point $P^* (x_1, \frac{1}{b}(\beta_2 x_1-r_2), \frac{m-cpx_1}{cpx_1})$ exists, if $\beta_2 x_1 > r_2$ and $m > cpx_1$, where $x_1$ is a root of the following polynomial

Here, we denote $A_i, i = 0, 1, 2, 3, 4$ such as

It is obvious to show that $A_4$ and $A_1$ has positive and negative signs, respectively, under the condition that we have

The number of positive real roots of (4.1) can be determined from the signs of $A_3$ and $A_2$. The number of roots will be revealed using the Descartes Rule. A number of positive roots are illustrated in Table 1.

Therefore, results can be determined by the following Lemma.

Lemma 4.1. The fractional prey-predator model $(1.2)$ has

1) A unique positive equilibrium point if cases 1–3 are satisfied.

2) More than one positive equilibrium when case 4 is satisfied.

Theorem 4.2. The trivial equilibrium point $P_0(0, 0, 0)$ is always unstable.

Proof. The characteristic equation of $J$ at the trivial equilibrium point $P_0(0, 0, 0)$ has the following form

The eigenvalues are $\lambda_1 = r_1-d; \lambda_2 = r_2$ and $\lambda_3 = -m$. Thus $|arg(\lambda_1)| = 0 < \frac{\alpha \pi}{2}$, whenever $r_1 > d$ and $|arg(\lambda_2)| = 0 < \frac{\alpha \pi}{2}$, since $r_2 > 0$ and $|arg(\lambda_3)| = \pi > \frac{\alpha \pi}{2}$. Hence, the convergence of the Mittag-Leffler function ensures that $P_0$ is always a saddle point. The instability of the trivial equilibrium point ensures that we do not expect a total extinction of the food chain cycle in the habitat.

Theorem 4.3. The prey equilibrium point $P_1(\frac{r_1-d}{a1}, 0, 0)$ is always unstable.

Proof. The characteristic equation of $J$ at the point $P_1$ has the following form

The eigenvalues are $\lambda_1 = r_1-d; \lambda_2 = r_2 -\frac{\beta_2(r_1 -d)}{a_1}$ and $\lambda_3 = \frac{cp (r_1 -d)}{a_1} -m$. Now, if $r_1 < d; r_2 a_1 < \beta_2(r_1 -d)$ and $cp(r_1 -d) < a_1m$, then $|arg(\lambda_1)| = \pi > \frac{\alpha \pi}{2}; |arg(\lambda_2)| = \pi > \frac{\alpha \pi}{2}$ and $|arg(\lambda_3)| = \pi > \frac{\alpha \pi}{2}$. Hence the prey equilibrium $P_1$ is locally asymptotically stable if $r_1 < d; r_2 a_1 < \beta_2(r_1 -d)$ and $cp(r_1 -d) < a_1 m$, which is not possible because the condition of existence for this point is $r_1 > d$. Thus, $P_1$ is unstable.

Theorem 4.4. The local asymptotic stability of $P_2(\bar{x_1}, \bar{x_2}, 0)$ holds if $C_1 > 0, C_3 > 0$ and $C_1 C_2 > C_3$.

Proof. The characteristic equation of the Jacobian matrix $J$ at $P_2(\bar{x_1}, \bar{x_2}, 0)$ is given by

where $C_1 = -J_{11} -J_{33} - J_{22}, C_2 = J_{22}J_{33}+ J_{11}J_{33}+J_{12}J_{21}+J_{13}J_{31}+J_{11}J_{22}, C_3 = -J_{11}J_{22}J_{33}-J_{12}J_{21}J_{33} -J_{13}J_{31}J_{22}, J_{11} = r_1 -(d + 2ax_1 + \beta_1x_2)$, $J_{12} = -\beta_1x_1, J_{13} = [ (1-\eta)/\kappa r1 + p]x_1, J_{21} = -\beta_2 x_2, J_{22} = r2-\beta_2x_1-a_2x_2, J_{31} = c(p+by)y$ and $J_{33} = cpx_1-m$. The above equation has three values $|arg(\lambda_1)| = \pi > \frac{\alpha \pi}{2}, |arg(\lambda_2)| = \pi > \frac{\alpha \pi}{2}, |arg(\lambda_3)| = \pi > \frac{\alpha \pi}{2}$, if $C_1 > 0, C_3 > 0$ and $C_1C_2 > C_3$ holds.

Theorem 4.5. The positive equilibrium point $P^*$ is conditionally locally asymptotically stable.

Proof. The characteristic equation associated with $J = [\rho_{ij}]_{i:j = 1, 2, 3}$ around the positive equilibrium point $P^*$ is

where

Let $D(\psi)$ be the discriminant of $\psi(\lambda)$, it can be given as

Then, we have

The positive equilibrium point $P^*$ is locally asymptotically stable provided one of the following satisfy:

1) $D(\psi) > 0, A_1 > 0$, $A_1 A_2 > A_3$ and $A_3 > 0$.

2) $D(\psi) < 0, A_1\ge 0, A_3 > 0, A_2\ge 0$ and $\theta < 2/3$.

3) $D(\psi) < 0, A_1 > 0, A_1 A_2 = A_3, A_2 > 0$ and for all $\theta \in (0, 1)$.

Theorem 4.6. The positive equilibrium point $P^*$ is conditionally globally asymptotically stable.

Proof. Let $V$ be the Lyapunov function defined such as

where $w_1$ and $w_2$ are non-negative constants. Taking fractional-order derivative on both sides, we have

Calculations show that we obtain

If we choose now $w_1 = \frac{\beta_1}{\beta_2}$ and $w_2 = \frac{1}{c}$, then it is seen that we have

Considering the ecological, environmental event, we assumed that all the parameters are positive and in addition, if $\eta < 1$, $x_1 > x_1^*$ and $y > y_1^*$, then we have $D^{\alpha} V(t)\leq 0$. This completes the proof.

5.   Bifurcation analysis

5.1.   Neimark-Sacker bifurcation

In this part, the Neimark-Sacker bifurcation conditions of the Caputo fractional order model are investigated. Now, let us take the following Caputo fractional-order model

We suppose that $E^{*}$ is an equilibrium point of system (1.2). Then (1.2) undergoes a Neimark-Sacker bifurcation around the point $E^{*}$ concerning the parameter $a$ at $a = a^{*}$ provided that the following conditions are satisfied:

$(i)$ the Jacobian matrix of system (1.2) at the equilibrium point $E^{*}$ has a pair of complex conjugate eigenvalues $\lambda_{1, 2} = \alpha_{j} \pm i \omega_{j}$ become purely imaginary at $a = a^{*}$.

$(ii) \; \xi_{1, 2}\left(\alpha, a^{*}\right) = 0$.

$(iii) \; \left.\frac{\partial \xi_{1, 2}}{\partial a}\right|_{a = a^{*}} \neq 0$.

where $\xi_{i}(\alpha, a) = -\min _{i = 1, 2}\left|\arg \left(\lambda_{i}(a)\right)\right|+\frac{\alpha \pi}{2}$.

The fractional-order derivative has an important role in increasing the stability of the considered model. Therefore, the conditions for which the system (1.2) undergoes Neimark-Sacker bifurcation concerning $\alpha$ is established as follows:

$(i)$ the Jacobian matrix of system (1.2) at $E^{*}$ has a pair of complex conjugate eigenvalues $\lambda_{1, 2} = \alpha_{j} \pm i \omega_{j}$ become purely imaginary at $\alpha = \alpha^{*}$.

$(ii) \; \phi_{1, 2}\left(\alpha^{*}\right) = 0$.

$(iii) \; \left.\frac{\partial \phi_{1, 2}}{\partial \alpha}\right|_{\alpha = \alpha^{*}} \neq 0$ where $\phi_{i}(\alpha) = \frac{\alpha \pi}{2}-\min _{i = 1, 2}\left|\arg \left(\lambda_{i}(\alpha)\right)\right|$.

5.2.   Neimark-Sacker bifurcation with discretization

The discretization of system $(1.2)$ is as follows

Then the solution of system $(1.2)$ for $t\in[0, h), \frac{t}{h}\in[0, 1)$ is given by

If we repeat the discretization process $n$ times, it is obtained that

Now for $t\in[nh, (n+1)h)$ and $t\rightarrow (n+1)h$ while $\alpha\rightarrow 1$, it is obtained that

The Jacobian matrix of system (5.1) around the interior equilibrium point $E^*$ is given by

where

While the characteristic equation of the Jacobian matrix (5.5) is given by

where

Now, using the Jury condition ^[9], the unique positive equilibrium $(x_1^*, x_2^*, y^*)$ is locally asymptotically stable if the following conditions are satisfied

To study the Neimark-Sacker bifurcation in the system (5.1), we need the following explicit criterion of Hopf bifurcation.

Lemma 5.3. (See ^[38]) Consider an $n$-dimensional discrete dynamical system $Y_{k+1} = f_\mu(Y_k)$ where $\mu\in R$ is the bifurcation parameter. Let $Y^*$ be a fixed point of $f_\mu$ and the characteristic polynomial for Jacobian matrix $J(Y^*) = (b_{ij})_{n\times n}$ of the $n$-dimensional map $f_\mu$ is given by

where $b_i = b_i(\mu, u), i = 1, 2, ..., n$ and $u$ is a control parameter to be determined. Let $\Delta_0^{\pm}(\mu, u) = 1, \Delta_1^{\pm}(\mu, u), ..., \Delta_n^{\pm}(\mu, u)$ be a sequence of determinants defined by $\Delta_i^{\pm}(\mu, u) = det(N_1\pm N_2), i = 1, 2, ..., n$ where

Moreover, the following conditions hold:

$(H1)$ Eigenvalue assignment: $\Delta_{n-1}^{-}(\mu_0, u) = 0, \Delta_{n-1}^{+}(\mu_0, u) > 0, Q_{\mu_0}(1) > 0, (-1)^nQ_{\mu_0}(-1) > 0, \Delta_i^{\pm}(\mu_0, u) > 0, i = n-3, n-5, ..., 1(or 2)$ when $n$ is even or odd respectively.

$(H2)$ Transversality condition: $[(\frac{d}{d\mu})(\Delta_{n-1}^{-}(\mu, u))]_{\mu = \mu_0}\neq 0$.

$(H3)$ Nonresonance condition: $cos(2\pi/j)\neq \phi$ or resonance condition $cos(2\pi/j) = \phi$ for $j = 3, 4, 5, ....$ and $\phi = -1+0.5Q_{\mu_0}(1)\Delta_{n-3}^{-}(\mu_0, u)/{\Delta^+_{n-2}(\mu_0, u)}$. Then Neimark-Sacker bifurcation happen if we take $\mu$ as a bifurcation parameter.

Theorem 5.4. The unique interior equilibrium point of model (5.4) undergoes Neimark-Sacker bifurcation if the following conditions hold:

where $B_1, B_2, B_3$ are described in Eq (5.6).

Proof. According to Lemma 1 for $n = 3$ we have in Eq (5.6), the characteristic polynomial of system (5.4) evaluated at this positive interior equilibrium. In this case, we obtain the following conclusions:

6.   Numerical simulations

In this part, we plot some simulations to support the previous results. We take $\alpha = 1.0, 0.9$ and $0.7$, to show the effect of the fractional-order derivatives on the dynamics of the proposed model. The system's initial conditions are chosen as $(0.14, 0.25, 0.37)$. The system parameter values are summarized in Table 2.

By analyzing the obtained figures, we reached some numerical fndings. In Figures 1–6, we have implemented the dynamics of the three species concerning time $t$. The behaviors of $x_1$, $x_2$ and $y$ populations are investigated with different parameter values listed in Table 2. We have also plotted the evolution of $x_1$, $x_2$ and $y$ concerning the fractional-order derivatives $\alpha = 1.0, 0.9$ and $0.7$. From Figure 1, we observe the stability of the prey equilibrium point for the parameter values of 2. 2 shows the existence of oscillations in the dynamics of the three species. In Figures 3 and 4, we take the fair parameter $k = 0.2$ and $k = 2$. It was obvious that the fair changed the behaviors of the species. In Figures 5 and 6, we choose the order of the fractional derivative $\alpha = 1$, $\alpha = 0.9$ and $\alpha = 0.7$. We noticed that the approach of the fractional-order derivative $\alpha$ to $0$ showed stability in the dynamical behavior of the three species.

7.   Conclusions

This work dealt with an ecological food chain cycle of a fractional-order predator-prey system. Both preys competed, while the single predator species showed cooperation in hunting species $x_1$. This collective strategy leads to a fear effect in the $x_1$ compartment. The second prey species $x_2$ was hunted individually so that an extreme change could not be noticed on this site. We used the Caputo fractional-order derivative considering the historical state of memory effect in the system. This study shows that fractional-order derivatives had a crucial role in controlling the stability of solutions in the three species. Furthermore, it affected the dynamics of solution, as can be seen in Figures 5 and 6. We observe from these figures that the approach to zero of the fractional-order derivative increases stability in the system of the three species $x_1$, $x_2$ and $y$. Moreover, the change of fair parameter $k$ showed infuence in temporal behaviors of the three species. In Figures 1 and 4, one can see that the fear of $x_1$ is effective when the population density of the predator increases. In this stage, the density of $x_1$ decreases. The case of predator-free equilibrium point shows that the compartment of $x_1$ increases the population's carrying capacity without any fear. However, $x_2$ is below the threshold of its compartment, and therefore it is extinct. $y$ doesn't have enough food in the habitat and shrinks.

Acknowledgments

The author T. Abdeljawad would like to thank Prince Sultan University for the support through the research lab TAS.

Conflict of interest

The authors declare no conflict of interest.