A Fractional-Order Eco-Epidemiological Leslie–Gower Model with Double Allee Effect and Disease in Predator

Abstract

In this paper, a fractional order of a modified Leslie–Gower predator-prey model with disease and the double Allee effect in predator population is proposed. Then, we analyze the important mathematical features of the proposed model such as the existence and uniqueness as well as the non-negativity and boundedness of solutions to the fractional-order system. Moreover, the local and global asymptotic stability conditions of all possible equilibrium points are investigated using Matignon’s condition and by constructing a suitable Lyapunov function, respectively. Finally, numerical simulations are presented to verify the theoretical results. We show numerically the occurrence of two limit cycles simultaneously driven by the order of the derivative, the bistability phenomenon for both the weak and strong Allee effect cases, and more dynamic behaviors such as the forward, backward, and saddle-node bifurcations which are driven by the transmission rate. We have found that the risk of extinction for the predator with a strong Allee effect is much higher when the spread of disease is relatively high.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

The study of dynamical behavior of the predator-prey model has been an important theme in mathematical biology; see, for instance [1–5], and references therein. The results of these studies could give us substantial elucidation about the qualitative behavior of the density of each population in the future time, such as stability, bifurcations, and chaos, without having any experimental laboratory. In the last few decades, special attention has been paid to the merger of the predator-prey model with the model for transmissible disease as a new branch of mathematical biology known as the eco-epidemiological model. It becomes an important tool in analyzing the effect of infectious diseases on the existence (extinction) of prey or/and predator populations. Several researchers have discussed the eco-epidemiological model incorporating ecological factors, such as harvesting [6], prey refuge [7, 8], social behavior [9], and the Alee effect [10].

One interesting phenomenon was introduced by [11] which demonstrated a condition where, at low densities of population, the presence of conspecifics could increase the per capita growth rate of the population. This phenomenon is called the Allee effect. Some population models involving the Allee effect have been developed to exhibit the important impact of this phenomenon on many aspects of population ecology, such as conservation of threatened species [12, 13], controlling pest species [14], and harvesting management [15, 16]. In general, there are two types of the Allee effects: the strong Allee effect and the weak Allee effect. The strong Allee effect has a threshold size population which is called the Allee threshold. When the density of population is under this threshold then the per capita growth rate of the population becomes negative. On the other hand, the weak Allee effect gives a reduction in the per capita growth rate when the population is low. Furthermore, several researchers investigate the dynamical behavior of a system incorporating any situation in which two or more components of the Allee effect can work simultaneously on a single population which is called the double Allee effect [1, 17]. These components could arrive from the reproduction or/and survival mechanisms such as the difficulty of finding a mate, cooperative breeding dependencies, cooperative antipredator behavior, and environmental conditioning, as shown in Table 1 [18].

Table 1

Parameter description in system (2).

Parameters	Ecological descriptions
${\widehat{r}}_1$	Intrinsic growth rate of prey
${\widehat{r}}_2$	Intrinsic growth rate of susceptible predator
$\widehat{K}$	Carrying capacity of prey
$\widehat{\theta }$	Capture rate of susceptible predator
$\widehat{b}$	Prey saturation constant
$\widehat{a}$	Another saturation constant
$\widehat{c}$	The strength of interference between susceptible predators in predation process
$\widehat{p}$	The Allee threshold
$\widehat{q}$	The auxiliary Allee effect constant
$\widehat{m}$	Intraspecific competition rate among predators
$\widehat{k}$	Alternative food for the predator
$\widehat{\beta }$	Transmission rate of the disease
$\widehat{\delta }$	Death rate of infected predator

Nowadays, fractional-order derivative modelling has become popular in many fields of science, such as physics, thermodynamics, biology, control theory, and many others (see, for example, [19–22]). Modeling using a fractional-order derivative has a major advantage in that it involves memory which comes from the fact that the fractional-order operators are nonlocal. More precisely, all of the previous conditions are captured in defining the fractional-order derivative. It has some advantages over the integer order derivative, especially the ability to describe the memory and hereditary properties which are inherent in various processes [23–25]. Some scholars had proven that systems with memory are more consistent and adequate with the real phenomena [26, 27].

Recently, [28] studied a fractional-order predator-prey model by assuming that the predation follows the Beddington-DeAngelis functional response and the double Allee effect on the predator population. By the last assumption, there is always a bistability condition for the strong Allee effect case. In other words, there is an extinction or existing condition for the predator depending on the choice of the initial condition. In this paper, we consider a predator-prey model in [28] by adding the assumption that there is a transmission of disease in the predator population. Biologically, the growth rates of a population must depend on the history of its previous conditions not only on the local conditions, and for this reason, the proposed model will use the fractional-order derivative which has a memory effect, to make it more accurate in predicting the future condition of the population. To the best of our knowledge, no investigation has been carried out on the dynamics of the proposed system. An important objective of this paper is to study how the spread of disease may affect changes in predator densities over time, where the intrinsic growth rate of the predator is affected by the double Allee effect.

We have arranged this paper in the following manner. First, we propose a model and give some useful preliminaries which are used in our analysis briefly. Next, we prove the basic properties of solutions. We also study the existence of the equilibrium points, their local stability, and the Hopf bifurcation. The global stability conditions for the equilibrium points are also examined. Moreover, some numerical simulations with a set of hypothetical parameters are demonstrated to validate the theoretical results and further explore the role of the disease in the predator-prey interaction. Finally, conclusions are given in the last section.

2. Formulation Model and Useful Preliminaries

In [28], the author’s considered a predator-prey model with memory effect, which is a modification of the Leslie–Gower model using the Caputo fractional derivative. The model abandoned the mass conservation principle held by the Lotka–Volterra model because of the assumption that both prey and predator obey the logistic law. They also assume that the predation process is governed by the Beddington-DeAngelis type functional response [29] which is a nonlinear prey-predator-dependent functional response, and that the intrinsic growth rate of the predator population is influenced by the double Allee effect, as given in the following form:$$\begin{array}{c}\text{(1)}& {}^C{{\displaystyle D}}_t^{\alpha }N={\widehat{r}}_{1}N1-\frac{N}{\widehat{K}}-\frac{\widehat{\theta }NP}{1+\widehat{c}N+\widehat{q}P},{}^C{{\displaystyle D}}_t^{\alpha }P=P\widehat{s}\frac{P-\widehat{p}}{P+\widehat{q}}-\frac{\widehat{m}P}{\widehat{k}+N}.\end{array}$$

In this paper, we assume that the predator population is split into two subpopulations, namely the susceptible predator $S$ and infected predator $I$ to include the disease transmission among predator. It is assumed that only susceptible predators can reproduce according to the logistic law, where the intrinsic growth rate of the susceptible predator is induced by the double Allee effect. The prey is only captured by the susceptible predator. This is reasonable from an ecological view as the infected predators are physically weak and lose their skill to catch food [30, 31]. It is also assumed that the susceptible predator has intraspecific competitive interactions with both the susceptible and infected predators. The transmission of the nonrecoverable disease in the predator population in the manner of the simple law of mass action and is vertically transferred from the susceptible predator to the infected predator compartment with a constant transmission rate $\beta$. Here, we also assume that the disease does not spread from the predator population to the prey population.

Considering the above assumptions, we construct a fractional-order model with the double Allee effect and disease in predator as follows:$$\begin{array}{c}\text{(2)}& {}^C{{\displaystyle D}}_t^{\alpha }N={\widehat{r}}_{1}N1-\frac{N}{\widehat{K}}-\frac{\widehat{\theta }NS}{\widehat{a}+\widehat{b}N+\widehat{c}S},{}^C{{\displaystyle D}}_t^{\alpha }DS=S{\widehat{r}}_2\frac{S-\widehat{p}}{S+\widehat{q}}-\frac{\widehat{m}S+\widehat{I}}{\widehat{k}+N}-\widehat{\beta }SI,\\ {}^C{{\displaystyle D}}_t^{\alpha }I=\widehat{\beta }SI-\widehat{\delta }I,\end{array}$$where ${}^C{{\displaystyle D}}_t^{\alpha }$ represents the Caputo fractional derivative with $0<\alpha <1$. The variables $N,S$, and $I$ denote the density of the prey, susceptible predator, and infected predator populations, respectively, and all parameters are shown in Table 1.

Since the right hand side of system (2) have time dimension ${\text{time}}^{-1}$ which is inconsistent with time dimension on the left hand side ${\text{time}}^{-\alpha }$, we modify the above system (2) as follows:$$\begin{array}{c}\text{(3)}& {}^C{{\displaystyle D}}_t^{\alpha }N={\widehat{r}}_1^{\alpha }N1-\frac{N}{\widehat{K}}-\frac{{\widehat{\theta }}^{\alpha }NS}{\widehat{a}+\widehat{b}N+\widehat{c}S},{}^C{{\displaystyle D}}_t^{\alpha }S=S{\widehat{r}}_2^{\alpha }\frac{S-\widehat{p}}{S+\widehat{q}}-\frac{{\widehat{m}}^{\alpha }S+I}{\widehat{k}+N}-{\widehat{\beta }}^{\alpha }SI,\\ {}^C{{\displaystyle D}}_t^{\alpha }I={\widehat{\beta }}^{\alpha }SI-{\widehat{\delta }}^{\alpha }I.\end{array}$$

For the sake of convenience, we redefine parameters of the system (3) as $r_1={\widehat{r}}_1^{\alpha },K=\widehat{K},\theta ={\widehat{\theta }}^{\alpha },a=\widehat{a},b=\widehat{b},c=\widehat{c},r_2={\widehat{r}}_2^{\alpha },p=\widehat{p},q=\widehat{q},m={\widehat{m}}^{\alpha },k=\widehat{k},\beta ={\widehat{\beta }}^{\alpha }$, and $\delta ={\widehat{\delta }}^{\alpha }$. Therefore, the system (3) with the initial conditions can be written in the following form:$$\begin{array}{c}\text{(4)}& {}^C{{\displaystyle D}}_t^{\alpha }N=r_{1}N1-\frac{N}{K}-\frac{\theta NS}{a+bN+cS},{}^C{{\displaystyle D}}_t^{\alpha }S=S{r}_2\frac{S-p}{S+q}-\frac{mS+I}{k+N}-\beta SI,\\ {}^C{{\displaystyle D}}_t^{\alpha }I=\beta SI-\delta I,\end{array}$$where $N0>0,S0>0$, and $I0>0$. Next, we briefly introduce definition and some basic properties of the Caputo fractional derivative.

Definition 1.

(see [32]). Let $\alpha \in \mathrm{0,1}$. The Caputo fractional derivative for the real valued function $ht$ is defined as follows:$$\begin{array}{}\text{(5)}& {}^C{{\displaystyle D}}_t^{\alpha }ht=\frac{1}{\Gamma 1-\alpha }{\displaystyle {\int }_0^t\frac{h\prime s}{{t-s}^{\alpha }}}\text{d}s,\end{array}$$where $\mathrm{\Gamma }\cdot$ is the Gamma function, $t\ge 0$.

To verify the non-negativity and uniform boundedness of the solutions of the system (4), we require the following lemma for the fractional derivative.

Lemma 2 (see [33]).

Suppose that $ht$ is continuous in $a,b$ and its derivatives ${}^C{{\displaystyle D}}_t^{\alpha }ht$ is continuous in $a,b$, for $0<\alpha <1$, then we have$$\begin{array}{c}\text{(6)}& ht=ha+\frac{1}{\Gamma \alpha }{}^C{{\displaystyle D}}_t^{\alpha }h\tau {t-a}^{\alpha },a\le \tau \le t,\forall t\in a,b.\end{array}$$

Lemma 3 (see [33]).

Suppose that $ht$ is continuous in $a,b$ and ${}^C{{\displaystyle D}}_t^{\alpha }ht$ is continuous in $a,b$ for $0<\alpha <1$. If ${}^C{{\displaystyle D}}_t^{\alpha }ht\ge 0,\forall t\in a,b$, then $ht$ is a nondecreasing function for all $t\in a,b$. If ${}^C{{\displaystyle D}}_t^{\alpha }ht\le 0,\forall t\in a,b$, then $ht$ is a nonincreasing function for all $t\in a,b$.

Lemma 4 (comparison lemma [34]).

Suppose $ht$ is continuous in $0,\infty$. If $ht$ satisfies$$\begin{array}{}\text{(7)}& {}^C{{\displaystyle D}}_t^{\alpha }ht\le -\phi ht+\mu ,\hspace{1em}h0=h_0\in ℝ,\end{array}$$where $0<\alpha <1,\phi ,\mu \in {ℝ}^2,\text{and}\phi \ne \text{0}$, then$$\begin{array}{}\text{(8)}& ht\le h_0-\frac{\mu }{\phi }{E}_{\alpha }-\phi t^{\alpha }+\frac{\mu }{\phi }.\end{array}$$

Next, we use the following theorem to determine the local stability behavior of the equilibrium points of the system (4).

Theorem 5 (Matignon condition [35]).

Consider the following Caputo fractional-order system with initial value$$\begin{array}{}\text{(9)}& {}^C{{\displaystyle D}}_t^{\alpha }\stackrel{⟶}{h}t=\stackrel{⟶}{f}t,\stackrel{⟶}{h},\hspace{1em}\stackrel{⟶}{h}0={\stackrel{⟶}{h}}_0,\end{array}$$where $\stackrel{⟶}{h}\in {ℝ}^n$ and $0<\alpha <1$. If ${\stackrel{⟶}{h}}^{\ast }$satisfies $\stackrel{⟶}{f}t,{\stackrel{⟶}{h}}^{\ast }=\stackrel{⟶}{0}$, then ${\stackrel{⟶}{h}}^{\ast }$ is an equilibrium point. If all eigenvalues ${\lambda }_i$, $i=1,\dots ,n$ of the Jacobian matrix $J{\stackrel{⟶}{h}}^{\ast }$ satisfy $\text{arg}{\lambda }_i>\alpha \pi /2$ then the equilibrium points ${\stackrel{⟶}{h}}^{\ast }$ is locally asymptotically stable.

To establish the conditions for the global stability of the equilibrium points of the system (4), we provide the following lemma.

Lemma 6 (Volterra-type Lyapunov function [36]).

Let $ht\in {ℝ}_{+}$ be a continuous and derivable function. Then, for $t\ge 0$,$$\begin{array}{c}\text{(10)}& {}^C{{\displaystyle D}}_t^{\alpha }ht-h^{\ast }-h^{\ast }\ln \frac{ht}{h^{\ast }}\le 1-\frac{h^{\ast }}{ht}{}^C{{\displaystyle D}}_t^{\alpha }ht,h^{\ast }\in {ℝ}_{+},\forall \alpha \in \mathrm{0,1}.\end{array}$$

Lemma 7 (generalized LaSalle invariance principle [37]).

Suppose $\mathrm{\Omega }$ is a bounded closed set. Every solution of ${}^C{{\displaystyle D}}_t^{\alpha }ht=fht$ which starts from a point in $\mathrm{\Omega }$ remains in $\mathrm{\Omega }$ for all $t\ge 0$. If there exists $Vh:\mathrm{\Omega }⟶ℝ$ with continuous first partial derivatives satisfying$$\begin{array}{}\text{(11)}& {}^C{{\displaystyle D}}_t^{\alpha }V{|}_{{}^C{{\displaystyle D}}_t^{\alpha }ht=fht}\le 0,\end{array}$$then every solution $ht$ starting in $\mathrm{\Omega }$ tends to $M$ as $t⟶\infty$ where $M$ is the largest invariant set of $E$ with$$\begin{array}{}\text{(12)}& E≔h{{}^C{{\displaystyle D}}_t^{\alpha }V}_{{}^C{{\displaystyle D}}_t^{\alpha }ht=fht}=0.\end{array}$$

3. Basic Properties of Solutions

In this section, we present some basic properties of solutions of the system (4) such as the non-negativity, boundedness, existence, and uniqueness of solutions.

3.1. Non-Negativity and Boundedness of Solutions

From the biological point of view, we only concern to system (4) in the non-negative and bounded solutions. We consider ${ℝ}_{+}^3≔N,S,I\in {ℝ}^3:N\ge 0,S\ge 0,I\ge 0$ as the set of all non-negative solutions. The following theorem shows the non-negativity and boundedness of the solutions of system (4).

Theorem 8.

All solutions of the system (4) which initiate in ${ℝ}_{+}^3$ are non-negative and uniformly bounded.

Proof.

Let $N0>0,S0>0$, and $I0>0$ and assume that $Nt\ge 0$ for $t\in 0,\infty$. Suppose the assumption is not true, then there exists $t^{\ast }$ which satisfies $0\le t<t^{\ast }$ such that $Nt>0$ for $t\in 0,t^{\ast }$, $N{t}^{\ast }=0$, and $Nt<0$ for $t>t^{\ast }$. From the first equation of system (4), we have$$\begin{array}{}\text{(13)}& {}^C{{\displaystyle D}}_t^{\alpha }{N{t}^{\ast }|}_{N{t}^{\ast }=0}=0.\end{array}$$

According to Lemma 3, we get $Nt=0$ for $t>t^{\ast }$ which contradicts to $Nt<0,\forall t>t^{\ast }$. Therefore, we have $Nt\ge 0$ for $t\in 0,\infty$. Using the similar way, it can be shown that $St\ge 0$ and $It\ge 0$ for $t\in 0,\infty$. Hence, all solutions of the system (4) are non-negative.

Next, we show the boundedness of solutions of the system (4) by using the fractional comparison lemma. From the first equation of system (4), we obtain$$\begin{array}{c}\text{(14)}& {}^C{{\displaystyle D}}_t^{\alpha }Nt+Nt=r_{1}N1-\frac{N}{K}-\frac{\theta NS}{a+bN+cS}+N\\ \le -r_1\frac{N^2}{K}+r_1+1N\\ =-\frac{r_1}{K}{N-\frac{K{r}_1+1}{2r_1}}^2\\ +\frac{K{r_1+1}^2}{4r_1}\\ \le \frac{K{r_1+1}^2}{4r_1}.\end{array}$$

Based on Lemma 4, we have$$\begin{array}{c}\text{(15)}& Nt\le N0-\frac{K{r_1+1}^2}{4r_1}{E}_{\alpha }-t^{\alpha }+\frac{K{r_1+1}^2}{4r_1},\end{array}$$where $E_{\alpha }$ is the Mittag–Leffler function. Lemma 5 and Corollary 6 in [38] implies that $E_{\alpha }-t^{\alpha }⟶0$ as $t⟶\infty$, then it follows$$\begin{array}{}\text{(16)}& 0\le Nt\le \frac{K{r_1+1}^2}{4r_1}={\eta }_1,\hspace{1em}\text{as\hspace{0.17em}}t⟶\infty .\end{array}$$

Thus, for any $\epsilon >0$, $Nt$ with initial condition in ${ℝ}_{+}^3$ is uniformly bounded to the region$$\begin{array}{}\text{(17)}& {\Theta }_1=Nt\le {\eta }_1+\epsilon ,\epsilon >0.\end{array}$$

Define $Wt=St+It$ and adding the last two equations of system (4). Then, we get$$\begin{array}{}\text{(18)}& {}^C{{\displaystyle D}}_t^{\alpha }Wt\le r_{2}S-p{r}_2\frac{S}{S+q}-m\frac{S^2}{k+N}-\delta I.\end{array}$$

For $\xi >0$ and $Nt\le {\eta }_1+\epsilon$, we have$$\begin{array}{c}\text{(19)}& {}^C{{\displaystyle D}}_t^{\alpha }Wt+\xi Wt\le -\frac{m}{k+{\eta }_1+\epsilon }{S-\frac{r_2+\xi k+{\eta }_1+\epsilon }{2m}}^2\\ +\frac{{r_2+\xi }^{2}k+{\eta }_1+\epsilon }{4m}\\ +\xi -\delta I-p{r}_2\frac{S}{S+q}.\end{array}$$

If we choose $\xi <\delta$, then$$\begin{array}{c}\text{(20)}& {}^C{{\displaystyle D}}_t^{\alpha }Wt+\xi Wt\le -\frac{m}{k+{\eta }_1+\epsilon }{S-\frac{r_2+\xi k+{\eta }_1+\epsilon }{2m}}^2+\frac{{r_2+\xi }^{2}k+{\eta }_1+\epsilon }{4m}-p{r}_2\frac{S}{S+q}.\\ \le \frac{{r_2+\xi }^{2}k+{\eta }_1+\epsilon }{4m}-p{r}_2\frac{S}{S+q}.\end{array}$$

Using the same argument as before and letting $\epsilon ⟶0$, then for the weak Alee effect $p<0$ we have the following equation:$$\begin{array}{c}\text{(21)}& 0\le Wt\le \frac{{r_2+\xi }^{2}k+{\eta }_1}{4\xi m}-\frac{p{r}_2}{\xi q}={\eta }_2,\hspace{1em}\text{as\hspace{0.17em}}t⟶\infty .\end{array}$$

On the other hand, for the strong Alee effect $p>0$ we have the following equation:$$\begin{array}{c}\text{(22)}& 0\le Wt\le \frac{{r_2+\xi }^{2}k+{\eta }_1}{4\xi m}={\eta }_3,\hspace{1em}\text{as\hspace{0.17em}}t⟶\infty .\end{array}$$

Thus, for any $\epsilon >0$, the combination of the solutions $St$ and $It$ with initial conditions started in ${ℝ}_{+}^3$ are uniformly bounded in the region$$\begin{array}{}\text{(23)}& {\Theta }_2=\begin{array}{c}St+It\le {\eta }_2+\epsilon ,\epsilon >0,p<0,\hfill \\ St+It\le {\eta }_3+\epsilon ,\epsilon >0,p>0.\hfill \end{array}\end{array}$$

3.2. Existence and Uniqueness of the Solution

To investigate the existence and uniqueness of solution of system (4), we apply the locally Lipschitz condition in the region $0,T\times \mathrm{\Upsilon }$, where $\mathrm{\Upsilon }=N,S,I\in {ℝ}^3:\text{max}N,S,I\le \mathrm{\Phi }$, for sufficient large $\mathrm{\Phi }$ and $T<\infty$.

Theorem 9.

For each $Z0=N0,S0,I0\in \mathrm{\Upsilon }$, there exists a unique solution $Zt=Nt,St,It\in \mathrm{\Upsilon }$ of system (4) with initial condition $Z0$, which is defined for all $t\ge 0$.

Proof.

For any $Z,\overline{Z}\in \mathrm{\Upsilon }$, we consider a mapping of a continuous function $FZ=F_{1}Z,F_{2}Z,F_{3}Z$ where$$\begin{array}{c}\text{(24)}& F_{1}Z=r_{1}N-\frac{r_1{N}^2}{K}-\frac{\theta NS}{a+bN+cS},F_{2}Z=\frac{r_2{S}^2}{S+q}-\frac{p{r}_{2}S}{S+q}-\frac{m{S}^2}{k+N}-\frac{mSI}{k+N}-\beta SI,\\ F_{3}Z=\beta SI-\delta I.\end{array}$$

It follows that$$\begin{array}{c}\text{(25)}& FZ-F\overline{Z}=F_{1}Z-F_1\overline{Z}+F_{2}Z-F_2\overline{Z}\\ +F_{3}Z-F_3\overline{Z}\\ =|r_{1}N-\frac{r_1{N}^2}{K}-\frac{\theta NS}{a+bN+cS}\\ -r_1\overline{N}-\frac{r_1{\overline{N}}^2}{K}-\frac{\theta \overline{N}\overline{S}}{a+b\overline{N}+c\overline{S}}|\\ +|\frac{r_2{S}^2}{S+q}-\frac{p{r}_{2}S}{S+q}-\frac{m{S}^2}{k+N}-\frac{mSI}{k+N}\\ -\beta SI-\frac{r_2{\overline{S}}^2}{\overline{S}+q}-\frac{p{r}_2\overline{S}}{\overline{S}+q}-\frac{m{\overline{S}}^2}{k+\overline{N}}\\ -\frac{m\overline{S}\overline{I}}{k+\overline{N}}-\beta \overline{S}\overline{I}|\\ +\beta SI-\delta I-\beta \overline{S}\overline{I}-\delta \overline{I}\\ \le r_{1}N-\overline{N}+\frac{r_1}{K}{N}^2-{\overline{N}}^2\\ +\theta \frac{NS}{a+bN+cS}-\frac{\overline{N}\overline{S}}{a+b\overline{N}+c\overline{S}}\\ +r_2\frac{S^2}{S+q}-\frac{{\overline{S}}^2}{\overline{S}+q}\\ +p{r}_2\frac{S}{S+q}-\frac{\overline{S}}{\overline{S}+q}\\ +m\frac{S^2}{k+N}-\frac{{\overline{S}}^2}{k+\overline{N}}\\ +m\frac{SI}{k+N}-\frac{\overline{S}\overline{I}}{k+\overline{N}}+\beta SI-\overline{S}\overline{I}\\ +\beta SI-\overline{S}\overline{I}+\delta I-\overline{I}\\ \le r_{1}N-\overline{N}+\frac{r_1}{K}N+\overline{N}N-\overline{N}\\ +\frac{\theta aN+b\overline{N}N}{a^2}S-\overline{S}\\ +\frac{\theta a\overline{S}+c\overline{S}S}{a^2}N-\overline{N}\\ +\frac{r_2}{q^2}\overline{S}SS-\overline{S}+\frac{r_2}{q}S+\overline{S}S-\overline{S}\\ +\frac{p{r}_2}{q}S-\overline{S}+\frac{mk+\overline{N}S+\overline{S}}{k^2}S-\overline{S}\\ +\frac{m{\overline{S}}^2}{k^2}N-\overline{N}+\frac{mk+\overline{N}I}{k^2}S-\overline{S}\\ +\frac{m\overline{N}S}{k^2}I-\overline{I}+\frac{m\overline{S}I}{k^2}N-\overline{N}\\ +2\beta SI-\overline{I}+2\beta IS-\overline{S}+\delta I-\overline{I}\\ \le M_{1}N-\overline{N}+M_{2}S-\overline{S}+M_{3}I-\overline{I}\\ \le MZ-\overline{Z},\end{array}$$where$$\begin{array}{c}\text{(26)}& M_1=r_1+\frac{2r_1}{K}+\frac{\theta }{a}\Phi +\frac{\theta c}{a^2}+\frac{2m}{k^2}{\Phi }^2,M_2=\frac{p{r}_2}{q}+\frac{\theta }{a}+\frac{2r_2}{q}+\frac{3m}{k}+2\beta \Phi \\ +\frac{\theta b}{a^2}+\frac{r_2}{q^2}+\frac{3m}{k^2}{\Phi }^2,\\ M_3=\delta +2\beta \Phi +\frac{m}{k^2}{\Phi }^2,\\ M=\max M_1,M_2,M_3.\end{array}$$

Since $FZ$ satisfies the Lipschitz condition with respect to $Z$, then system (4) with initial condition $Z0$ has a unique solution $Zt\in \mathrm{\Upsilon }$ for all $t\ge 0$. The proof is completed.

3.3. Existence and Stability of Equilibrium Points

To study the dynamical behaviors of system (4), we first investigate the equilibrium points of the system (4) which are the constant solutions to the following system:$$\begin{array}{c}\text{(27)}& r_{1}N1-\frac{N}{K}-\frac{\theta NS}{a+bN+cS}=0,S{r}_2\frac{S-p}{S+q}-\frac{mS+I}{k+N}-\beta SI=0,\\ \beta SI-\delta I=0.\end{array}$$

Next, the local stability of equilibrium points of system (4) is evaluated by computing the eigenvalues of the Jacobian matrix of system (4) at the equilibrium point $N^{\ast },S^{\ast },I^{\ast }$ which is given as follows:$$\begin{array}{}\text{(28)}& J{N}^{\ast },S^{\ast },I^{\ast }=\begin{array}{ccc}{J}_{\mathrm{1,1}}& J_{\mathrm{1,2}}& 0\\ J_{\mathrm{2,1}}& J_{\mathrm{2,2}}& J_{\mathrm{2,3}}\\ 0& J_{\mathrm{3,2}}& J_{\mathrm{3,3}}\end{array},\end{array}$$where$$\begin{array}{c}\text{(29)}& J_{\mathrm{1,1}}=r_{1}1-\frac{N^{\ast }}{K}-\frac{r_1{N}^{\ast }}{K}-\frac{\theta S^{\ast }}{a+b{N}^{\ast }+c{S}^{\ast }}+\frac{b\theta N^{\ast }{S}^{\ast }}{{a+b{N}^{\ast }+c{S}^{\ast }}^2},\\ J_{\mathrm{1,2}}=-\frac{\theta N^{\ast }}{a+b{N}^{\ast }+c{S}^{\ast }}+\frac{c\theta N^{\ast }{S}^{\ast }}{{a+b{N}^{\ast }+c{S}^{\ast }}^2}\\ J_{\mathrm{2,1}}=\frac{m{S}^{\ast }{S}^{\ast }+I^{\ast }}{{k+N^{\ast }}^2}\\ J_{\mathrm{2,2}}=r_2\frac{S^{\ast }-p}{S^{\ast }+q}-\frac{m{S}^{\ast }+I^{\ast }}{k+N^{\ast }}+\\ S^{\ast }\frac{r_2}{S^{\ast }+q}-\frac{r_2{S}^{\ast }-p}{{S^{\ast }+q}^2}-\frac{m}{k+N^{\ast }}-\beta I^{\ast }\\ J_{\mathrm{2,3}}=-\frac{m{S}^{\ast }}{k+N^{\ast }}+\beta S^{\ast }\\ J_{\mathrm{3,2}}=\beta I^{\ast }\\ J_{\mathrm{3,3}}=\beta S^{\ast }-\delta .\end{array}$$

3.4. Equilibrium Analysis in the Case of the Weak Allee Effect in Predator

In this subsection, the existence and stability of the equilibrium points of the system (4) for the weak Alee effect $p<0$ are given. By solving system (9) simultaneously, all positive equilibrium points of system (4) are obtained as follows:

(a) Trivial equilibrium point: $E_0=\mathrm{0,0,0}$, which always exists.

(b) Axial equilibrium point: $E_N=K,\mathrm{0,0}$, which always exists; and $E_S=0,S_w,0$, where $S_w=k{r}_2-mq/2m+\sqrt{{k{r}_2-mq/2m}^2-kp{r}_2/m}$, which always exists.

(c) Planar equilibrium point: $E_{SI}=0,{\overline{S}}_w,{\overline{I}}_w$ where ${\overline{S}}_w=\delta /\beta$ and ${\overline{I}}_w=1/m+\beta k\beta k{r}_2\delta -\beta p-\delta m\delta +\beta q/\beta \delta +\beta q$ exists if $p<\text{min}0,\delta /\beta 1-m\delta +\beta q/\beta k{r}_2$; and $E_{NS}={\tilde{N}}_w,{\tilde{S}}_w,0$ exists if $p>{\tilde{S}}_{w}1-m{\tilde{S}}_w+q/k{r}_2$, where ${\tilde{N}}_w=m{\tilde{S}}_w{\tilde{S}}_w+q/r_2{\tilde{S}}_w-p-k$ and ${\tilde{S}}_w$ are all positive roots of the following quartic equation:$$\begin{array}{}\text{(30)}& l_1{S}^4+l_2{S}^3+l_3{S}^2+l_{4}S+l_5=0,\end{array}$$

  where$$\begin{array}{c}\text{(31)}& l_1=mb+c{r}_{2}m{r}_1,l_2=-bk-cp+a{r}_2+mbqm{r}_1\\ -mb+c{r}_{2}k+K{r}_2-mq{r}_1+\theta K{r}_2^2,\\ l_3=-p-bk+am{r}_1{r}_2--bk-cp+a{r}_2\\ +mbqk+K{r}_2-mq{r}_1\\ +mb+c{r}_{2}pk+K{r}_2{r}_1-2\theta K{r}_2^{2}p,\\ l_4=p{r}_{2}k+K{r}_2-mq-bk+a{r}_1\\ +-bk-cp+a{r}_2+mbq\\ pk+K{r}_2{r}_1+\theta K{r}_2^2{p}^2,\\ l_5=-p-bk+apk+K{r}_2{r}_2{r}_1.\end{array}$$

(d) Interior equilibrium point: $E^{\ast }=N_w^{\ast },S_w^{\ast },I_w^{\ast }$, where

$$\begin{array}{c}\text{(32)}& N_w^{\ast }=\frac{bK-a+c{S}_w^{\ast }}{2b}\pm \sqrt{{\Delta }_w},S_w^{\ast }=\frac{\delta }{\beta },I_w^{\ast }=\frac{r_2{S}_w^{\ast }-pk+N_w^{\ast }-m{S}_w^{\ast }{S}_w^{\ast }+q}{S_w^{\ast }+qm+\beta k+N_w^{\ast }},\\ {\Delta }_w={\frac{bK-a+c{S}_w^{\ast }}{2b}}^2+\frac{K{r}_{1}a+c{S}_w^{\ast }-\theta S_w^{\ast }}{b{r}_1}.\end{array}$$

Assume that $k+N_w^{\ast }>\delta m\delta +\beta q/\beta r_2\delta -\beta p$, then we have the following results:

(a) If ${\mathrm{\Delta }}_w<0$, then system (4) has no interior equilibrium point.

(b) If ${\mathrm{\Delta }}_w=0$ and $a\beta +c\delta <bK\beta$, then system (4) has a unique interior equilibrium point $E^{\ast }=N_w^{\ast },\delta /\beta ,I_w^{\ast }$ where $N_w^{\ast }=bK\beta -a\beta +c\delta /2b\beta$ and $I_w^{\ast }=\beta r_2\delta -\beta pk+N_w^{\ast }-\delta m\delta +\beta q/\beta \delta +\beta qm+\beta k+N_w^{\ast }$.

(c) If ${\mathrm{\Delta }}_w>0$ and $a\beta +c\delta <\text{min}bK\beta ,\delta \theta /r_1$, then system (4) has two positive interior equilibrium points $E_{\mathrm{1,2}}^{\ast }=N_{wi}^{\ast },\delta /\beta ,I_{wi}^{\ast }$, $i=\mathrm{1,2}$, where

$$\begin{array}{c}\text{(33)}& N_{w_{\mathrm{1,2}}}^{\ast }=\frac{bK\beta -a\beta +c\delta }{2b\beta }\pm \sqrt{{\Delta }_w},I_{w_{\mathrm{1,2}}}^{\ast }=\frac{\beta r_2\delta -\beta pk+{\widehat{N}}_{w_{\mathrm{1,2}}}-\delta m\delta +\beta q}{\beta \delta +\beta qm+\beta k+{\widehat{N}}_{w_{\mathrm{1,2}}}},\\ {\Delta }_w={\frac{bK\beta -a\beta +c\delta }{2b\beta }}^2\\ +\frac{K{r}_{1}a\beta +c\delta -\delta \theta }{b{r}_1\beta }.\end{array}$$

Theorem 10.

The local stability of trivial, axial, and planar equilibrium points of system (4) for the weak Alee effect $p<0$ are summarized as follows.

(a) $E_0=\mathrm{0,0,0}$ is always a saddle point.

(b) $E_N=K,\mathrm{0,0}$ is always a saddle point.

(c) $E_S=0,S_w,0$ is locally asymptotically stable if $r_1<\theta S_w/a+c{S}_w$, $p+q<m/k{r}_2{S_w+q}^2$, and $S_w<\delta /\beta$.

(d) Let ${\widehat{\mathrm{\Delta }}}_w={\beta \delta r_{2}p+q/{\delta +\beta q}^2-m\delta /\beta k}^2-4\delta /\beta k\beta k{r}_2\delta -\beta p-m\delta \delta +\beta q/\delta +\beta q$. $E_{SI}=0,{\overline{S}}_w,{\overline{I}}_w$ is locally asymptotically stable if $r_1<\delta \theta /a\beta +c\delta$ and one of the following conditions holds:

(i) $\beta r_{2}p+q/{\delta +\beta q}^2<m/\beta k<r_2\delta -\beta p/\delta \delta +\beta q$, or

(ii) $\beta r_{2}p+q/{\delta +\beta q}^2>m/\beta k$, ${\widehat{\mathrm{\Delta }}}_w<0$, and $\alpha <{\alpha }^{\ast }=2/\pi {\tan }^{-1}\sqrt{{\widehat{\mathrm{\Delta }}}_w}/\beta r_{2}p+q/{\delta +\beta q}^2-m/\beta k$.

(e) Suppose that$$\begin{array}{c}\text{(34)}& {\omega }_1=-\frac{r_1{\tilde{N}}_w}{K}1-\frac{b{r}_{1}K-{\tilde{N}}_w}{\theta K{\tilde{S}}_w}+\frac{r_2{{\tilde{S}}_w-p}^2-pp+q}{{{\tilde{S}}_w+q}^2},{\omega }_2=\frac{r_1{\tilde{N}}_w}{K}1-\frac{b{r}_{1}K-{\tilde{N}}_w}{\theta K{\tilde{S}}_w}\\ \frac{r_2{{\tilde{S}}_w-p}^2-pp+q}{{{\tilde{S}}_w+q}^2}+\frac{\theta {\tilde{N}}_{w}a+b{\tilde{N}}_w}{{a+b{\tilde{N}}_w+c{\tilde{S}}_w}^2}\frac{m{\tilde{S}}_w^2}{{k+{\tilde{N}}_w}^2}.\end{array}$$

  $E_{NS}={\tilde{N}}_w,{\tilde{S}}_w,0$ is locally asymptotically stable if ${\tilde{S}}_w<\delta /\beta$ and

(i) ${\omega }_1<0$ and ${\omega }_2>0$, or

(ii) ${\omega }_1>0$, ${\omega }_1^2-4{\omega }_2<0$ and $\alpha <{\alpha }^{\ast }=2/\pi {\tan }^{-1}\sqrt{{\omega }_1^2-4{\omega }_2}/{\omega }_1$

Proof.

(a) In view of (28), around $E_0=\mathrm{0,0,0}$ we have$$\begin{array}{}\text{(35)}& J{E}_0=\begin{array}{ccc}{r}_1& 0& 0\\ 0& -\frac{r_{2}p}{q}& 0\\ 0& 0& -\delta \end{array}.\end{array}$$

  The corresponding eigenvalues of (35) are ${\lambda }_1=r_1>0,{\lambda }_2=-r_{2}p/q>0$, and ${\lambda }_3=-\delta <0$. Thus, $\text{arg}{\lambda }_{\mathrm{1,2}}=0<\alpha \pi /2$ and $\text{arg}{\lambda }_3=\pi >\alpha \pi /2$. Therefore, $E_0$ is always a saddle point.

(b) The Jacobian matrix in (28) around $E_N=K,\mathrm{0,0}$ is$$\begin{array}{}\text{(36)}& J{E}_N=\begin{array}{ccc}-r_1& -\frac{\theta K}{a+bK}& 0\\ 0& -\frac{r_{2}p}{q}& 0\\ 0& 0& -\delta \end{array}.\end{array}$$

  Then, the corresponding eigenvalues of (36) are ${\lambda }_1=-r_1<0,{\lambda }_2=-r_{2}p/q>0$, and ${\lambda }_3=-\delta <0$. Since $\text{arg}{\lambda }_2=0<\alpha \pi /2$ and $\text{arg}{\lambda }_{\mathrm{1,3}}=\pi >\alpha \pi /2$, therefore, $E_N$ is a saddle point.

(c) By substituting $E_S=0,S_w,0$ to the Jacobian matrix in (28), we have$$\begin{array}{}\text{(37)}& J{E}_S=\begin{array}{ccc}{a}_{11}& 0& 0\\ a_{21}& a_{22}& a_{23}\\ 0& 0& a_{33}\end{array},\end{array}$$

  where$$\begin{array}{c}\text{(38)}& a_{11}=r_1-\frac{\theta S_w}{a+c{S}_w},a_{21}=\frac{m{S}_w^2}{k^2},\\ a_{22}=\frac{k{r}_2{S}_{w}p+q-m{S}_w{S_w+q}^2}{k{S_w+q}^2},\\ a_{23}=-S_w\frac{m}{k}+\beta ,\\ a_{33}=\beta S_w-\delta .\end{array}$$

  The corresponding eigenvalues of (37) are ${\lambda }_1=r_1-\theta S_w/a+c{S}_w$, ${\lambda }_2=k{r}_2{S}_{w}p+q-m{S}_w{S_w+q}^2/k{S_w+q}^2$, and ${\lambda }_3=\beta S_w-\delta$. Therefore, $\text{arg}{\lambda }_{\mathrm{1,2,3}}=\pi >\alpha \pi /2$, whenever $r_1<\theta S_w/a+c{S}_w$, $p+q<m/k{r}_2{S_w+q}^2$, and $\delta >\beta S_w$.

(d) The Jacobian matrix in (28) around $E_{SI}=0,{\overline{S}}_w,{\overline{I}}_w$ is as follows:$$\begin{array}{}\text{(39)}& J{E}_{SI}=\begin{array}{ccc}{b}_{11}& 0& 0\\ b_{21}& b_{22}& b_{23}\\ 0& b_{32}& 0\end{array},\end{array}$$

  where$$\begin{array}{c}\text{(40)}& b_{11}=r_1-\frac{\delta \theta }{a\beta +c\delta },b_{21}=\frac{m\delta }{\beta k}\frac{{\delta }^2+\beta q+r_2\delta -\beta r_{2}p}{m+\beta k\delta +\beta q},\\ b_{22}=\frac{\beta \delta r_{2}p+q}{{\delta +\beta q}^2}-\frac{m\delta }{\beta k},\\ b_{23}=-\delta \frac{m}{\beta k}+1,\\ b_{32}=\frac{\beta k{r}_2\delta -\beta p-m\delta \delta +\beta q}{m+\beta k\delta +\beta q}.\end{array}$$

  The corresponding eigenvalues of (39) are ${\lambda }_1=r_1-\delta \theta /a\beta +c\delta$ and roots of the following equation:$$\begin{array}{c}\text{(41)}& {\lambda }^2-\frac{\beta \delta r_{2}p+q}{{\delta +\beta q}^2}-\frac{m\delta }{\beta k}\lambda +\frac{\delta }{\beta k}\frac{\beta k{r}_2\delta -\beta p-m\delta \delta +\beta q}{\delta +\beta q}=0.\end{array}$$

  Equation (41) has the following eigenvalues:$$\begin{array}{}\text{(42)}& {\lambda }_{\mathrm{2,3}}=\frac{1}{2}\frac{\beta \delta r_{2}p+q}{{\delta +\beta q}^2}-\frac{m\delta }{\beta k}\pm \sqrt{{\widehat{\Delta }}_w},\end{array}$$

  where$$\begin{array}{c}\text{(43)}& {\widehat{\Delta }}_w={\frac{\beta \delta r_{2}p+q}{{\delta +\beta q}^2}-\frac{m\delta }{\beta k}}^2-4\frac{\delta }{\beta k}\frac{\beta k{r}_2\delta -\beta p-m\delta \delta +\beta q}{\delta +\beta q}.\end{array}$$

  Notice that if $r_1<\delta \theta /a\beta +c\delta$, then $\text{arg}{\lambda }_1=\pi >\alpha \pi /2$. Therefore, the stability of $E_{SI}$ depends on ${\lambda }_{\mathrm{2,3}}$. If $\beta r_{2}p+q/{\delta +\beta q}^2<m/\beta k<r_2\delta -\beta p/\delta \delta +\beta q$, then $\text{arg}{\lambda }_{\mathrm{2,3}}>\alpha \pi /2$. Furthermore, if $\beta r_{2}p+q/{\delta +\beta q}^2>m/\beta k$ and ${\widehat{\mathrm{\Delta }}}_w<0$, then ${\lambda }_{\mathrm{2,3}}$ is a pair of complex conjugate eigenvalues with positive real parts. Thus, $\text{arg}{\lambda }_{\mathrm{2,3}}>\alpha \pi /2$ is achieved only if $\alpha <{\alpha }^{\ast }=2/\pi {\tan }^{-1}\sqrt{{\widehat{\mathrm{\Delta }}}_w}/\beta r_{2}p+q/{\delta +\beta q}^2-m/\beta k$. Using the Matignon’s condition (see Theorem 5), the theorem is completely proven.

(e) The Jacobian matrix (10) calculated at $E_{NS}={\tilde{N}}_w,{\tilde{S}}_w,0$ is given by

$$\begin{array}{}\text{(44)}& J{E}_{NS}=\begin{array}{ccc}{c}_{11}& c_{12}& 0\\ c_{21}& c_{22}& c_{23}\\ 0& 0& c_{33}\end{array},\end{array}$$where$$\begin{array}{c}\text{(45)}& c_{11}=-\frac{r_1{\tilde{N}}_w}{K}1-\frac{b{r}_{1}K-{\tilde{N}}_w}{\theta K{\tilde{S}}_w},c_{12}=-\frac{\theta {\tilde{N}}_{w}a+b{\tilde{N}}_w}{{a+b{\tilde{N}}_w+c{\tilde{S}}_w}^2},\\ c_{21}=\frac{m{\tilde{S}}_w^2}{{k+{\tilde{N}}_w}^2},\\ c_{22}=-\frac{r_2{{\tilde{S}}_w-p}^2-pp+q}{{{\tilde{S}}_w+q}^2},\\ c_{23}=-{\tilde{S}}_w\frac{m}{k+{\tilde{N}}_w}+\beta ,\\ c_{33}=\beta {\tilde{S}}_w-\delta .\end{array}$$

The eigenvalues of (44) are ${\lambda }_1=\beta {\tilde{S}}_w-\delta$ and the roots of the quadratic equation ${\lambda }^2-{\omega }_1\lambda +{\omega }_2=0$ where$$\begin{array}{c}\text{(46)}& {\omega }_1=-\frac{r_1{\tilde{N}}_w}{K}1-\frac{b{r}_{1}K-{\tilde{N}}_w}{\theta K{\tilde{S}}_w}+\frac{r_2{{\tilde{S}}_w-p}^2-pp+q}{{{\tilde{S}}_w+q}^2},\\ {\omega }_2=\frac{r_1{\tilde{N}}_w}{K}1-\frac{b{r}_{1}K-{\tilde{N}}_w}{\theta K{\tilde{S}}_w}\\ \frac{r_2{{\tilde{S}}_w-p}^2-pp+q}{{{\tilde{S}}_w+q}^2}\\ +\frac{\theta {\tilde{N}}_{w}a+b{\tilde{N}}_w}{{a+b{\tilde{N}}_w+c{\tilde{S}}_w}^2}\frac{m{\tilde{S}}_w^2}{{k+{\tilde{N}}_w}^2}.\end{array}$$

If ${\tilde{S}}_w<\delta /\beta$, then $\text{arg}{\lambda }_1=\pi >\alpha \pi /2$. Therefore, the stability of $E_{NS}$ depends on ${\lambda }_{\mathrm{2,3}}$. If ${\omega }_1<0$ and ${\omega }_2>0$ then $\text{arg}{\lambda }_{\mathrm{2,3}}>\alpha \pi /2$. Furthermore, if ${\omega }_1>0$ and ${\omega }_1^2<4{\omega }_2$, then ${\lambda }_{\mathrm{2,3}}$ is a pair of complex conjugate eigenvalues. Thus, $\text{arg}{\lambda }_{\mathrm{2,3}}>\alpha \pi /2$ is attained if $\alpha <{\alpha }^{\ast }=2/\pi {\tan }^{-1}\sqrt{{\omega }_1^2-4{\omega }_2}/{\omega }_1$. So, by the Matignon’s condition (see Theorem 5) system (4) exhibits locally asymptotically stable behavior around $E_{NS}$.

Theorem 11.

Stability condition of interior equilibrium point for weak Alee effect $p<0$. Suppose that$$\begin{array}{c}\text{(47)}& {\eta }_1=\frac{r_1{N}_w^{\ast }}{K}1-\frac{b{r}_1{K-N_w^{\ast }}^2}{\theta K{S}_w^{\ast }}-S_w^{\ast }\frac{r_{2}p+q}{{S_w^{\ast }+q}^2}-\frac{m}{k+N_w^{\ast }},\\ {\eta }_2=-S_w^{\ast }\frac{r_1{N}_w^{\ast }}{K}1-\frac{b{r}_1{K-N_w^{\ast }}^2}{\theta K{S}_w^{\ast }}\\ \frac{r_{2}p+q}{{S_w^{\ast }+q}^2}-\frac{m}{k+N_w^{\ast }}\\ +\beta I_w^{\ast }\frac{m}{k+N_w^{\ast }}+\beta \\ -\frac{m\theta N_w^{\ast }a+b{N}_w^{\ast }{S}_w^{\ast }+I_w^{\ast }}{{a+b{N}_w^{\ast }+c{S}_w^{\ast }}^2{k+N_w^{\ast }}^2},\\ {\eta }_3=\frac{\beta r_1{N}_w^{\ast }{S}_w^{\ast }{I}_w^{\ast }}{K}1-\frac{b{r}_1{K-N_w^{\ast }}^2}{\theta K{S}_w^{\ast }}\\ \frac{m}{k+N_w^{\ast }}+\beta ,\\ D_{w}P=18{\eta }_1{\eta }_2{\eta }_3+{{\eta }_1{\eta }_2}^2-4{\eta }_3{\eta }_1^3\\ -4{\eta }_2^3-27{\eta }_3^2.\end{array}$$

$E^{\ast }=N_w^{\ast },S_w^{\ast },I_w^{\ast }$ is locally asymptotically stable if

(i) $D_{w}P>0,{\eta }_1>0,{\eta }_3>0$, and${\eta }_1{\eta }_2>{\eta }_3$.

(ii) $D_{w}P<0,{\eta }_1\ge 0,{\eta }_2\ge 0,{\eta }_3>0$, and$0<\alpha <2/3$.

(iii) $D_{w}P<0,{\eta }_1<0,{\eta }_2<0$, and$2/3<\alpha <1$.

(iv) $D_{w}P<0,{\eta }_1>0,{\eta }_2>0$,${\eta }_1{\eta }_2={\eta }_3$, and$0<\alpha <1$.

Proof.

The Jacobian matrix (10) evaluated at interior equilibrium point $E^{\ast }=N_w^{\ast },S_w^{\ast },I_w^{\ast }$ is given by$$\begin{array}{}\text{(48)}& J{E}^{\ast }=\begin{array}{ccc}{d}_{11}& d_{12}& 0\\ d_{21}& d_{22}& d_{23}\\ 0& d_{32}& 0\end{array},\end{array}$$where$$\begin{array}{c}\text{(49)}& d_{11}=-\frac{r_1{N}_w^{\ast }}{K}1-\frac{b{r}_1{K-N_w^{\ast }}^2}{\theta K{S}_w^{\ast }},d_{12}=-\frac{\theta N_w^{\ast }a+b{N}_w^{\ast }}{{a+b{N}_w^{\ast }+c{S}_w^{\ast }}^2},\\ d_{21}=\frac{m{S}_w^{\ast }{S}_w^{\ast }+I_w^{\ast }}{{k+N_w^{\ast }}^2},\\ d_{22}=S_w^{\ast }\frac{r_{2}p+q}{{S_w^{\ast }+q}^2}-\frac{m}{k+N_w^{\ast }},\\ d_{23}=-S_w^{\ast }\frac{m}{k+N_w^{\ast }}+\beta \\ d_{32}=\beta I_w^{\ast }.\end{array}$$

The corresponding eigenvalues of (48) are the roots of the cubic equation $P\lambda ={\lambda }^3+{\eta }_1{\lambda }^2+{\eta }_2\lambda +{\eta }_3=0$ where the discriminant $D_{w}P$ of the cubic equation $P\lambda$ is as follows:$$\begin{array}{}\text{(50)}& D_{w}P=18{\eta }_1{\eta }_2{\eta }_3+{{\eta }_1{\eta }_2}^2-4{\eta }_3{\eta }_1^3-4{\eta }_2^3-27{\eta }_3^2.\end{array}$$

If $D_{w}P>0$, all the roots of the cubic equation are real; and there is only one real root and two complex conjugate roots if $D_{w}P<0$. Therefore, the Routh-Hurwitz criterion for Caputo fractional-order [39] completes the proof of the stability condition for $E^{\ast }$.

3.5. Equilibrium Analysis in the Strong Allee Effect in Predator

In this subsection, the existence and stability of the equilibrium points of system (4) for the strong Alee effect $p>0$ are given. All positive equilibrium points of system (4) are as follows.

(a) Trivial equilibrium point: ${\mathrm{\Xi }}_0=\mathrm{0,0,0}$, which always exists,

(b) Axial equilibrium point: ${\mathrm{\Xi }}_N=K,\mathrm{0,0}$, which always exists; and ${\mathrm{\Xi }}_S=0,S_s,0$where $S_{s_{\mathrm{1,2}}}=k{r}_2-mq/2m\pm \sqrt{{k{r}_2-mq/2m}^2-kp{r}_2/m}$. Assume that $k{r}_2>mq$, the existence of ${\mathrm{\Xi }}_S$ is described as follows.

(i) If ${k{r}_2-mq/2m}^2<kp{r}_2/m$, then ${\mathrm{\Xi }}_S$ does not exist.

(ii) If ${k{r}_2-mq/2m}^2<kp{r}_2/m$, then there exists a unique ${\mathrm{\Xi }}_S=0,k{r}_2-mq/2m,0$

(iii) If ${k{r}_2-mq/2m}^2<kp{r}_2/m$, then there exists two${\mathrm{\Xi }}_S=0,k{r}_2-mq/2m\pm \sqrt{{k{r}_2-mq/2m}^2-kp{r}_2/m},0$

(c) Planar equilibrium point: ${\mathrm{\Xi }}_{SI}=0,{\overline{S}}_s,{\overline{I}}_s$ exists if $p<{\overline{S}}_{s}1-m\delta +\beta q/\beta k{r}_2$, where$$\begin{array}{c}\text{(51)}& {\overline{S}}_s=\frac{\delta }{\beta },\text{and}{\overline{I}}_w=\frac{1}{m+\beta k}\frac{\beta k{r}_2\delta -\beta p-\delta m\delta +\beta q}{\beta \delta +\beta q};\end{array}$$

  and ${\mathrm{\Xi }}_{NS}={\tilde{N}}_s,{\tilde{S}}_s,0$ exists if ${\tilde{S}}_{s}1-m{\tilde{S}}_s+q/k{r}_2<p<{\tilde{S}}_s$, where ${\tilde{N}}_s=m{\tilde{S}}_s{\tilde{S}}_s+q/r_2{\tilde{S}}_s-p-k$ and ${\tilde{S}}_s$ are also all positive roots of quartic (30).

(d) Interior equilibrium point: ${\mathrm{\Xi }}^{\ast }=N_s^{\ast },S_s^{\ast },I_s^{\ast }$, where

$$\begin{array}{c}\text{(52)}& N_{s_{\mathrm{1,2}}}^{\ast }=\frac{bK-a+c{S}_s^{\ast }}{2b}\pm \sqrt{{\Delta }_s},S_s^{\ast }=\frac{\delta }{\beta },\\ I_s^{\ast }=\frac{r_2{S}_s^{\ast }-pk+N_s^{\ast }-m{S}_s^{\ast }{S}_s^{\ast }+q}{S_s^{\ast }+qm+\beta k+N_s^{\ast }},\\ {\Delta }_s={\frac{bK-a+c{S}_s^{\ast }}{2b}}^2\\ +\frac{K{r}_{1}a+c{S}_s^{\ast }-\theta S_s^{\ast }}{b{r}_1}.\end{array}$$

Assume that $S_s^{\ast }>p$ and $k+N_s^{\ast }>m{S}_s^{\ast }{S}_s^{\ast }+q/r_2{S}_s^{\ast }-p$, then we have the following results.

(a) If ${\mathrm{\Delta }}_s<0$, then system (4) has no interior equilibrium point.

(b) If ${\mathrm{\Delta }}_s=0$ and $a\beta +c\delta <bK\beta$, then system (4) has a unique interior equilibrium point ${\mathrm{\Xi }}_1^{\ast }=N_s^{\ast },\delta /\beta ,I_s^{\ast }$ where $N_s^{\ast }=bK\beta -a\beta +c\delta /2b\beta$ and $I_s^{\ast }=\beta r_2\delta -\beta pk+N_s^{\ast }-\delta m\delta +\beta q/\beta \delta +\beta qm+\beta k+N_s^{\ast }$.

(c) If ${\mathrm{\Delta }}_s>0$ and $a\beta +c\delta <\text{min}bK\beta ,\delta \theta /r_1$, then system (4) has two positive interior equilibrium points ${\mathrm{\Xi }}_{\mathrm{1,2}}^{\ast }=N_{si}^{\ast },\delta /\beta ,I_{si}^{\ast }$ for $i=\mathrm{1,2}$, where

$$\begin{array}{c}\text{(53)}& N_{s_{\mathrm{1,2}}}^{\ast }=\frac{bK\beta -a\beta +c\delta }{2b\beta }\pm \sqrt{{\Delta }_s},I_{s_{\mathrm{1,2}}}^{\ast }=\frac{\beta r_2\delta -\beta pk+N_{s_{\mathrm{1,2}}}^{\ast }-\delta m\delta +\beta q}{\beta \delta +\beta qm+\beta k{N}_{s_{\mathrm{1,2}}}^{\ast }},\\ {\Delta }_s={bK\beta -a\beta +c\delta /2b\beta }^2\\ +K{r}_{1}a\beta +c\delta -\delta \theta /b{r}_1\beta .\end{array}$$

Theorem 12.

The local stability of all equilibrium points of system (4) for the strong Alee effect $p>0$ is summarized as follows.

(a) ${\mathrm{\Xi }}_0=\mathrm{0,0,0}$ is a saddle point.

(b) ${\mathrm{\Xi }}_N=K,\mathrm{0,0}$ is always locally asymptotically stable.

(c) ${\mathrm{\Xi }}_S=0,S_s,0$ is locally asymptotically stable if $r_1<\theta S_s/a+c{S}_s$, $p+q<m/k{r}_2{S_s+q}^2$, and $S_s<\delta /\beta$.

(d) Let ${\overline{S}}_s>p$ and ${\widehat{\mathrm{\Delta }}}_s={\beta \delta r_{2}p+q/{\delta +\beta q}^2-m\delta /\beta k}^2-4\delta /\beta k\beta k{r}_2\delta -\beta p-m\delta \delta +\beta q/\delta +\beta q$. ${\mathrm{\Xi }}_{SI}=0,{\overline{S}}_s,{\overline{I}}_s$ is locally asymptotically stable if $r_1<\delta \theta /a\beta +c\delta$ and

(i) $\beta r_{2}p+q/{\delta +\beta q}^2<m/\beta k<r_2\delta -\beta p/\delta \delta +\beta q$, or

(ii) $\beta r_{2}p+q/{\delta +\beta q}^2>m/\beta k$, ${\widehat{\mathrm{\Delta }}}_s<0$, and $\alpha <{\alpha }^{\ast }=2/\pi {\tan }^{-1}\sqrt{{\widehat{\mathrm{\Delta }}}_s}/\beta r_{2}p+q/{\delta +\beta q}^2-m/\beta k$.

(e) Suppose that$$\begin{array}{c}\text{(54)}& {\omega }_{s_1}=-\frac{r_1{\tilde{N}}_s}{K}1-\frac{b{r}_{1}K-{\tilde{N}}_s}{\theta K{\tilde{S}}_s}+\frac{r_2{{\tilde{S}}_s-p}^2-pp+q}{{{\tilde{S}}_s+q}^2},\\ {\omega }_{s_2}=\frac{r_1{\tilde{N}}_s}{K}1-\frac{b{r}_{1}K-{\tilde{N}}_s}{\theta K{\tilde{S}}_s}\\ \frac{r_2{{\tilde{S}}_s-p}^2-pp+q}{{{\tilde{S}}_s+q}^2}\\ +\frac{\theta {\tilde{N}}_{s}a+b{\tilde{N}}_s}{{a+b{\tilde{N}}_s+c{\tilde{S}}_s}^2}\frac{m{\tilde{S}}_s^2}{{k+{\tilde{N}}_s}^2}.\end{array}$$

  ${\mathrm{\Xi }}_{NS}={\tilde{N}}_s,{\tilde{S}}_s,0$ is locally asymptotically stable if ${\tilde{S}}_s<\delta /\beta$ and

(i) ${\omega }_{s_1}<0$, or

(ii) ${\omega }_{s_2}>0$,${\omega }_{s_1}^2-4{\omega }_{s_2}<0$, and$\alpha <{\alpha }^{\ast }=2/\pi {\tan }^{-1}\sqrt{{\omega }_{s_1}^2-4{\omega }_{s_2}}/{\omega }_{s_1}$.

(f) Suppose that

$$\begin{array}{c}\text{(55)}& {\xi }_1=\frac{r_1{N}_s^{\ast }}{K}1-\frac{b{r}_1{K-N_s^{\ast }}^2}{\theta K{S}_s^{\ast }}-S_s^{\ast }\frac{r_{2}p+q}{{S_s^{\ast }+q}^2}-\frac{m}{k+N_s^{\ast }},\\ {\xi }_2=-S_s^{\ast }\frac{r_1{N}_s^{\ast }}{K}1-\frac{b{r}_1{K-N_s^{\ast }}^2}{\theta K{S}_s^{\ast }}\\ \frac{r_{2}p+q}{{S_s^{\ast }+q}^2}-\frac{m}{k+N_s^{\ast }}\\ +\beta I_s^{\ast }\frac{m}{k+N_s^{\ast }}+\beta \\ -\frac{m\theta N_s^{\ast }a+b{N}_s^{\ast }{S}_s^{\ast }+I_s^{\ast }}{{a+b{N}_s^{\ast }+c{S}_s^{\ast }}^2{k+N_s^{\ast }}^2}\\ {\xi }_3=\frac{\beta r_1{N}_s^{\ast }{S}_s^{\ast }{I}_s^{\ast }}{K}1-\frac{b{r}_1{K-N_s^{\ast }}^2}{\theta K{S}_s^{\ast }}\\ \frac{m}{k+N_s^{\ast }}+\beta \\ D_{s}P=18{\xi }_1{\xi }_2{\xi }_3+{{\xi }_1{\xi }_2}^2-4{\xi }_3{{\xi }_1}^3\\ -4{{\xi }_2}^3-27{{\xi }_3}^2.\end{array}$$

${\mathrm{\Xi }}^{\ast }=N_s^{\ast },S_s^{\ast },I_s^{\ast }$ is locally asymptotically stable if

(i) $D_{s}P>0,{\xi }_1>0,{\xi }_3>0,$and${\xi }_1{\xi }_2>{\xi }_3$.

(ii) $D_{s}P<0,{\xi }_1\ge 0,{\xi }_2\ge 0,{\xi }_3>0$, and$0<\alpha <2/3$.

(iii) $D_{s}P<0,{\xi }_1<0,{\xi }_2<0$, and$2/3<\alpha <1$.

(iv) $D_{s}P<0,{\xi }_1>0,{\xi }_2>0$,${\xi }_1{\xi }_2={\xi }_3$, and$0<\alpha <1$.

Theorem 12 has similar proof to Theorem 10 and Theorem 11.

3.6. Hopf Bifurcation

In this subsection, we study the conditions of a Hopf bifurcation around the equilibrium point of system (4) when a parameter is varied. This bifurcation ensures a stability change when system (4) passes the critical value which coincides with the emergence of the limit cycle. Hopf bifurcation can occur both in the first-order systems and fractional-order systems. The fundamental difference between the two is the convergence of the limit set of solutions, known as the limit cycle, to the solution of the system. In the first-order systems, the limit cycle converges to the periodic solution; in the fractional-order systems, instead of converging to the periodic solution, the limit cycle converges to the periodic signal [40, 41].

Let us consider the following three-dimensional fractional order system:$$\begin{array}{}\text{(56)}& {}^C{{\displaystyle D}}_t^{\alpha }x=gx,\hspace{1em}0<\alpha \le 1,x\in {ℝ}^3.\end{array}$$

According to Theorem 5, the stability of the system dynamics is significantly affected by the order of the derivative $\alpha$. Suppose that $E_{\ast }$ is an equilibrium point of system (19) and ${\alpha }^{\ast }$ is the critical value of bifurcation parameter. The conditions under which system (19) undergoes a Hopf bifurcation driven by $\alpha$ near the equilibrium point $E_{\ast }$ are given in [42].

(1) The Jacobian matrix at $E_{\ast }$ has one real negative eigen value ${\lambda }_1<0$ and a pair of complex conjugate eigen values ${\lambda }_{\mathrm{2,3}}=\psi \pm i\varphi$ (where $\psi >0$),

(2) $\chi {\alpha }^{\ast }={\alpha }^{\ast }\pi /2-\text{arg}{\lambda }_{\mathrm{2,3}}\alpha =0$, and

(3) The transversality condition: ${\partial \chi \alpha /\partial \alpha |}_{\alpha ={\alpha }^{\ast }}\ne 0$.

Note that the critical value of ${\alpha }^{\ast }$ is the solution of $\chi {\alpha }^{\ast }={\alpha }^{\ast }\pi /2-\text{arg}{\lambda }_{\mathrm{2,3}}\alpha =0$, i.e., ${\alpha }^{\ast }=2/\pi {\tan }^{-1}\varphi /\psi$. Therefore, we have the following theorem.

Theorem 13.

Existence of Hopf bifurcation driven by $\alpha$ around $E_{SI}$, $E_{NS}$, or $E^{\ast }$ for the weak Allee effect case.

(i) Suppose that $r_1<\delta \theta /a\beta +c\delta$, $\beta r_{2}p+q/{\delta +\beta q}^2>m/\beta k$, and ${\tilde{\mathrm{\Delta }}}_w<0$. System (4) undergoes a Hopf bifurcation around $E_{SI}$ when $\alpha$ passes through$$\begin{array}{}\text{(57)}& {\alpha }^{\ast }=\frac{2}{\pi }{\tan }^{-1}\frac{\sqrt{{\tilde{\Delta }}_w}}{\beta r_{2}p+q/{\delta +\beta q}^2-m/\beta k}.\end{array}$$

(ii) If $\tilde{S}<\delta /\beta$, ${\omega }_1>0$, and ${\omega }_1^2-4{\omega }_2<0$. Then, system (4) undergoes a Hopf bifurcation around $E_{NS}$ when parameter bifurcation $\alpha$ passes through$$\begin{array}{}\text{(58)}& {\alpha }^{\ast }=\frac{2}{\pi }{\tan }^{-1}\frac{\sqrt{{\omega }_1^2-4{\omega }_2}}{{\omega }_1}.\end{array}$$

(iii) Suppose that characteristic equation $P\lambda ={\lambda }^3+{\eta }_1{\lambda }^2+{\eta }_2\lambda +{\eta }_3=0$ has one real negative eigen value and a pair of complex conjugate eigen values ${\lambda }_{\mathrm{2,3}}={\zeta }_1\pm i{\zeta }_2$ where ${\zeta }_1>0$. System (4) undergoes a Hopf bifurcation around the interior equilibrium point $E^{\ast }=N_w^{\ast },S_w^{\ast },I_w^{\ast }$ when $\alpha$ passes through

$$\begin{array}{}\text{(59)}& {\alpha }^{\ast }=2\pi {\tan }^{-1}\frac{\sqrt{{\zeta }_1^2-4{\zeta }_2}}{{\zeta }_1}.\end{array}$$

Remark 14.

In addition to the existence of Hopf bifurcation driven by order of the fractional derivative $\alpha$, the transmission rate $\beta$ also can be considered as a parameter bifurcation. Theorem 4 in [42] gives the existence condition for Hopf bifurcation by varying $\beta$. Since it is difficult to find out the exact critical values of $\beta ={\beta }^{\ast }$, therefore we numerically investigated the existence of Hopf bifurcation in system (4) driven by parameter $\beta$.

4. Global Stability

By using the appropriate Lyapunov function, we investigate the global stability of the stable equilibrium points of system (4), both for the weak Allee effect and the strong Allee effect cases in the predator.

4.1. For the Weak Allee Effect in Predator

Theorem 15.

$E_S=0,S_w,0$ is globally asymptotically stable, if$$\begin{array}{c}\text{(60)}& -q\frac{S_w}{{\eta }_2+q}-1\le p\le -\frac{{\eta }_2+q}{r_2{S}_w}\frac{r_{1}K}{4}+m{S}_w^2,S_w<\frac{\delta }{\beta },\end{array}$$hold.

Proof.

Let $Nt,St,It$ be any positive solution of system (4). From (20), an $\epsilon >0$ can be chosen such that$$\begin{array}{c}\text{(61)}& \frac{p{r}_2{S}_w}{{\eta }_2+\epsilon +q}\le -\frac{r_{1}K}{4}+m{S}_w^2,\frac{p}{q}\ge 1-\frac{S_w}{{\eta }_2+\epsilon +q},\\ S_w<\frac{\delta }{m+\beta }<\frac{\delta }{\beta }.\end{array}$$

Further, from Theorem 8, we have $0\le St\le {\eta }_2+\epsilon$. We now define a positive definite Lyapunov function at $E_S=0,S_w,0$ as follows:$$\begin{array}{}\text{(62)}& L_{1}N,S,I=N+S-S_w-S_w\text{ln}\frac{S}{S_w}+I.\end{array}$$

By using Lemma 6, the fractional time derivative of $L_{1}N,S,I$ along the solutions of system (4) is given by$$\begin{array}{c}\text{(63)}& {}^C{{\displaystyle D}}_t^{\alpha }{L}_{1}N,S,I\le {}^C{{\displaystyle D}}_t^{\alpha }N+\frac{S-S_w}{S}{}^C{{\displaystyle D}}_t^{\alpha }S+{}^C{{\displaystyle D}}_t^{\alpha }I\\ =r_{1}1-\frac{N}{K}-\frac{\theta S}{a+bN+cS}N\\ +S-S_w\frac{r_{2}S-p}{S+q}-\frac{mS+I}{k+N}-\beta I\\ +\beta SI-\delta I\\ =-\frac{r_1}{K}{N-\frac{K}{2}}^2+\frac{r_{1}K}{4}-\frac{\theta NS}{a+bN+cS}\\ +\frac{r_2{S}^2}{S+q}-\frac{p{r}_{2}S}{S+q}-\frac{r_2{S}_{w}S}{S+q}+\frac{p{r}_2{S}_w}{S+q}\\ -\frac{m{S}_{w}S}{k+N}+\frac{m{S}_w^2}{k+N}-\frac{m}{k+N}{S-S_w}^2\\ -\frac{mSI}{k+N}+\frac{m{S}_{w}I}{k+N}+\beta S_{w}I-\delta I.\\ \le \frac{r_{1}K}{4}+\frac{p{r}_2{S}_w}{{\eta }_2+\epsilon +q}+m{S}_w^2\\ +r_{2}1-\frac{p}{q}-\frac{S_w}{{\eta }_2+\epsilon +q}S\\ -\frac{m}{k+N}{S-S_w}^2+m{S}_w+\beta S_w-\delta I.\end{array}$$

If condition in (61) is satisfied then we obtain the following:$$\begin{array}{}\text{(64)}& {}^C{{\displaystyle D}}_t^{\alpha }{L}_{1}N,S,I\le -\frac{m}{k+N}{S-S_w}^2.\end{array}$$

In this case, ${}^C{{\displaystyle D}}_t^{\alpha }{L}_{1}N,S,I\le 0,\hspace{1em}\forall N,S,I\in {ℝ}_{+}^3$, and ${}^C{{\displaystyle D}}_t^{\alpha }{L}_{1}N,S,I=0$ at $E_S$. Based on that, the only invariant set on which ${}^C{{\displaystyle D}}_t^{\alpha }{L}_{1}N,S,I=0$ is the singleton $E_S$. Using Lemma 7, the sufficient condition for the global asymptotic stability of $E_S$ is achieved.

Theorem 16.

The planar equilibrium point $E_{SI}=0,{\overline{S}}_w,{\overline{I}}_w$ is globally asymptotically stable, if$$\begin{array}{c}\text{(65)}& -\frac{\beta {\eta }_2+q}{\delta p}\frac{r_{1}K}{4}+\frac{{\delta }^{2}m}{{\beta }^{2}k}+\delta {\overline{I}}_w+\frac{\delta m{\eta }_2}{\beta k}\le r_2\le \frac{p{r}_2}{q}+\frac{\delta r_2}{\beta {\eta }_2+q}+\frac{\delta m}{\beta {\eta }_1+k}+\beta {\overline{I}}_w,\end{array}$$holds.

Proof.

Assume that $Nt,St,It$ be the positive solution of system (4), so we can choose $\epsilon >0$ such that$$\begin{array}{c}\text{(66)}& \frac{\delta p{r}_2}{\beta {\eta }_2+\epsilon +q}\le -\frac{r_{1}K}{4}+\frac{{\delta }^{2}m}{{\beta }^{2}k}+\delta {\overline{I}}_w+\frac{\delta m}{\beta k}{\eta }_2+\epsilon \\ r_2\le \frac{p{r}_2}{q}+\frac{\delta r_2}{\beta {\eta }_2+\epsilon +q}\\ +\frac{\delta m}{\beta k+{\eta }_1+\epsilon }+\beta {\overline{I}}_w.\end{array}$$

From Theorem 8, we also have that $Nt$ and $St$ are bounded, i.e., $0\le Nt\le {\eta }_1+\epsilon$ and $0\le St\le {\eta }_2+\epsilon$. Now, we define a positive definite Lyapunov function at $E_{SI}=0,{\overline{S}}_w,{\overline{I}}_w$ as$$\begin{array}{c}\text{(67)}& L_{2}N,S,I=N+S-\frac{\delta }{\beta }-\frac{\delta }{\beta }\text{ln}\frac{\beta S}{\delta }+I-{\overline{I}}_w-{\overline{I}}_w\text{ln}\frac{I}{{\overline{I}}_w}.\end{array}$$

As before, we have the fractional time derivative of $L_{2}N,S,I$ along the solutions of system (4) as follows:$$\begin{array}{c}\text{(68)}& {}^C{{\displaystyle D}}_t^{\alpha }{L}_{2}N,S,I\le {}^C{{\displaystyle D}}_t^{\alpha }N+\frac{S-\delta /\beta }{S}{}^C{{\displaystyle D}}_t^{\alpha }S+\frac{I-{\overline{I}}_w}{I}{}^C{{\displaystyle D}}_t^{\alpha }I\\ =r_{1}1-\frac{N}{K}-\frac{\theta S}{a+bN+cS}N\\ +S-\frac{\delta }{\beta }\frac{r_{2}S-p}{S+q}-\frac{mS+I}{k+N}-\beta I\\ \le \frac{r_{1}K}{4}+\frac{r_2{S}^2}{S+q}-\frac{p{r}_{2}S}{S+q}-\frac{r_2\delta /\beta S}{S+q}+\frac{p{r}_2\delta /\beta }{S+q}\\ -\frac{m\delta /\beta S}{k+N}+\frac{m{\delta /\beta }^2}{k+N}-\frac{m{S-\delta /\beta }^2}{k+N}-\frac{mSI}{k+N}\\ +\frac{m\delta /\beta I}{k+N}-\beta {\overline{I}}_{w}S+{\overline{I}}_w\delta \\ \le \frac{r_{1}K}{4}+\frac{\delta p{r}_2}{\beta {\eta }_2+\epsilon +q}+\frac{{\delta }^{2}m}{{\beta }^{2}k}+\delta {\overline{I}}_w\\ +\frac{\delta m}{\beta k}{\eta }_2+\epsilon +r_2-\frac{p{r}_2}{q}\\ -\frac{\delta r_2}{\beta {\eta }_2+\epsilon +q}-\frac{\delta m}{\beta k+{\eta }_1+\epsilon }-\beta {\overline{I}}_{w}S\\ -\frac{m}{k+N}{S-\frac{\delta }{\beta }}^2-\frac{m}{k+N}SI.\end{array}$$

Observe that when the condition of (66) is fulfilled, then we have the following:$$\begin{array}{}\text{(69)}& {}^C{{\displaystyle D}}_t^{\alpha }{L}_{2}N,S,I\le -\frac{m}{k+N}{S-\frac{\delta }{\beta }}^2-\frac{m}{k+N}SI.\end{array}$$

One can easily show that ${}^C{{\displaystyle D}}_t^{\alpha }{L}_{2}N,S,I\le 0,\hspace{1em}\forall N,S,I\in {ℝ}_{+}^3$, and ${}^C{{\displaystyle D}}_t^{\alpha }{L}_{2}N,S,I=0$ at $E_{SI}$. Therefore, the only invariant set on which ${}^C{{\displaystyle D}}_t^{\alpha }{L}_{2}N,S,I=0$ is the singleton $E_{SI}$. Following Lemma 7, whenever $E_{SI}$ exists and$$\begin{array}{c}\text{(70)}& -\frac{\beta {\eta }_2+q}{\delta p}\frac{r_{1}K}{4}+\frac{{\delta }^{2}m}{{\beta }^{2}k}+\delta {\overline{I}}_w+\frac{\delta m{\eta }_2}{\beta k}\le r_2\le \frac{p{r}_2}{q}+\frac{\delta r_2}{\beta {\eta }_2+q}+\frac{\delta m}{\beta {\eta }_1+k}+\beta {\overline{I}}_w,\end{array}$$then it is globally asymptotically stable.

Theorem 17.

The planar equilibrium point of system (4) $E_{NS}={\tilde{N}}_w,{\tilde{S}}_w,0$ is globally asymptotically stable if$$\begin{array}{c}\text{(71)}& {\tilde{N}}_w\ge K,r_1\le \frac{{\tilde{S}}_w-p{r}_2/{\eta }_2+q-m{\tilde{S}}_w/k}{{\tilde{N}}_w{\tilde{N}}_w/K-1},r_2\le \frac{m{\tilde{S}}_w/k+{\eta }_1-\theta {\tilde{N}}_w/a}{1-{\tilde{S}}_w/{\eta }_2+q-p/q},\delta \ge \frac{m{\tilde{S}}_w}{k}+\beta {\tilde{S}}_w.\end{array}$$

Proof.

Let $Nt,St,It$ be any positive solution of system (4). Then, from (24), an $\epsilon >0$ can be chosen such that$$\begin{array}{c}\text{(72)}& 1-\frac{{\tilde{N}}_w}{K}<0,\frac{r_1{\tilde{N}}_w^2}{K}-r_1{\tilde{N}}_w<-\frac{p{r}_2{\tilde{S}}_w}{{\eta }_2+\epsilon +q}+\frac{m{\tilde{S}}_w^2}{k},\\ \frac{m{\tilde{S}}_w}{k}+\beta {\tilde{S}}_w<\delta \\ \frac{\theta {\tilde{N}}_w}{a}+r_2<\frac{r_2{\tilde{S}}_w}{{\eta }_2+\epsilon +q}+\frac{p{r}_2}{q}\\ +\frac{m{\tilde{S}}_w}{k+{\eta }_1+\epsilon }.\end{array}$$

Also, from Theorem 8, we have that $Nt$ and $St$ are bounded such that $0\le Nt\le {\eta }_1+\epsilon$ and $0\le St\le {\eta }_2+\epsilon$. Let us consider a positive definite Lyapunov function at $E_{NS}={\tilde{N}}_w,{\tilde{S}}_w,0$ as follows:$$\begin{array}{c}\text{(73)}& L_{3}N,S,I=N-{\tilde{N}}_w-{\tilde{N}}_w\text{ln}\frac{N}{{\tilde{N}}_w}+S-{\tilde{S}}_w-{\tilde{S}}_w\text{ln}\frac{N}{{\tilde{S}}_w}+I.\end{array}$$

Based on Lemma 6 and by calculating the fractional time derivative of $L_{3}N,S,I$ along the solutions of system (4),$$\begin{array}{c}\text{(74)}& {}^C{{\displaystyle D}}_t^{\alpha }{L}_{3}N,S,I\le \frac{N-{\tilde{N}}_w}{N}{}^C{{\displaystyle D}}_t^{\alpha }N+\frac{S-{\tilde{S}}_w}{S}{}^C{{\displaystyle D}}_t^{\alpha }S+{}^C{{\displaystyle D}}_t^{\alpha }I\\ =\frac{N-{\tilde{N}}_w}{N}{r}_{1}1-\frac{N}{K}-\frac{\theta S}{a+bN+cS}N\\ +\frac{S-{\tilde{S}}_w}{S}\frac{r_{2}S-p}{S+q}-\frac{mS+I}{k+N}-\beta IS\\ +\beta SI-\delta I\\ =-\frac{r_1}{K}{N-{\tilde{N}}_w}^2+r_{1}N-\frac{r_1{\tilde{N}}_w}{K}N-r_1{\tilde{N}}_w\\ +\frac{r_1{\tilde{N}}_w^2}{K}-\frac{\theta NS}{a+bN+cS}+\frac{\theta {\tilde{N}}_{w}S}{a+bN+cS}\\ +\frac{r_2{S}^2}{S+q}-\frac{r_{2}S{\tilde{S}}_w}{S+q}-\frac{p{r}_{2}S}{S+q}+\frac{p{r}_2{\tilde{S}}_w}{S+q}\\ -\frac{m}{k+N}{S-{\tilde{S}}_w}^2-\frac{m{\tilde{S}}_w}{k+N}S+\frac{m{\tilde{S}}_w^2}{k+N}\\ -\frac{mIS}{k+N}+\frac{mI{\tilde{S}}_w}{k+N}+\beta {\tilde{S}}_{w}I-\delta I\\ \le -\frac{r_1}{K}{N-{\tilde{N}}_w}^2+r_{1}1-\frac{{\tilde{N}}_w}{K}N\\ +\frac{r_1{\tilde{N}}_w^2}{K}-r_1{\tilde{N}}_w+\frac{p{r}_2{\tilde{S}}_w}{{\eta }_2+\epsilon +q}+\frac{m{\tilde{S}}_w^2}{k}\\ +\frac{\theta {\tilde{N}}_w}{a}+r_2-\frac{r_2{\tilde{S}}_w}{{\eta }_2+\epsilon +q}-\frac{p{r}_2}{q}\\ -\frac{m{\tilde{S}}_w}{k+{\eta }_1+\epsilon }S-\frac{m}{k+N}{S-{\tilde{S}}_w}^2\\ -\frac{m}{k+N}SI+\frac{m{\tilde{S}}_w}{k}+\beta {\tilde{S}}_w-\delta I.\end{array}$$

If condition in (72) is achieved, then we have the following:$$\begin{array}{c}\text{(75)}& {}^C{{\displaystyle D}}_t^{\alpha }{L}_{3}N,S,I\le -\frac{r_1}{K}{N-{\tilde{N}}_w}^2-\frac{m}{k+N}{S-{\tilde{S}}_w}^2-\frac{m}{k+N}SI.\end{array}$$

It is obvious that ${}^C{{\displaystyle D}}_t^{\alpha }{L}_{3}N,S,I\le 0,\hspace{1em}\forall N,S,I\in {ℝ}_{+}^3$, and ${}^C{{\displaystyle D}}_t^{\alpha }{L}_{3}N,S,I=0$ at $E_{NS}$. Hence, the singleton $E_{NS}$ is the only invariant set on which ${}^C{{\displaystyle D}}_t^{\alpha }{L}_{3}N,S,I=0$. Thus, we can conclude that $E_{NS}$ is globally asymptotically stable which allowed by Lemma 7. □

Theorem 18.

Suppose that$$\begin{array}{c}\text{(76)}& {\Omega }_1=\frac{2a{r}_1}{K}\frac{a\beta +b{N}_w^{\ast }\beta -c\delta }{a\beta +b{N}_w^{\ast }\beta +2b\delta }{\Omega }_2=\frac{2am}{k+{\eta }_{1}a+b{N}_w^{\ast }}\\ {\Omega }_3=\frac{\beta I_w^{\ast }q}{r_2}-q\\ {\Omega }_4=\frac{\delta m+\beta m{\eta }_2+{\beta }^{2}k{I}_w^{\ast }{\eta }_2+q}{\beta k{r}_2}.\end{array}$$

If $\theta \le \text{min}{\mathrm{\Omega }}_1,{\mathrm{\Omega }}_2$ and ${\mathrm{\Omega }}_3\le -p\le {\mathrm{\Omega }}_4$, then $E^{\ast }=N_w^{\ast },S_w^{\ast },I_w^{\ast }$ is globally asymptotically stable.

Proof.

Let $Nt,St,It$ be any positive solution of system (4), where $Nt$ and $St$ are bounded by Theorem 8. Further, one can use the following positive definite Lyapunov function to study the global stability of $E^{\ast }=N_w^{\ast },S_w^{\ast },I_w^{\ast }$:$$\begin{array}{c}\text{(77)}& L_{5}N,S,I=a+b{N}_w^{\ast }+c\frac{\delta }{\beta }N-N_w^{\ast }-N_w^{\ast }\ln \frac{N}{N_w^{\ast }}\\ +S-\frac{\delta }{\beta }-\frac{\delta }{\beta }\ln \frac{N}{\delta /\beta }+I-I_w^{\ast }-I_w^{\ast }\ln \frac{I}{I_w^{\ast }}.\end{array}$$

Again, using Lemma 6, we have the following equation:$$\begin{array}{c}\text{(78)}& {}^C{{\displaystyle D}}_t^{\alpha }{L}_{5}N,S,I\le a+b{N}_w^{\ast }+c\frac{\delta }{\beta }\frac{N-N_w^{\ast }}{N}{}^C{{\displaystyle D}}_t^{\alpha }N\\ +\frac{S-\delta /\beta }{S}{}^C{{\displaystyle D}}_t^{\alpha }S+\frac{I-I_w^{\ast }}{I}{}^C{{\displaystyle D}}_t^{\alpha }I\\ =a+b{N}_w^{\ast }+c\frac{\delta }{\beta }N-N_w^{\ast }\\ r_{1}1-\frac{N}{K}-\frac{\theta S}{a+bN+cS}\\ +S-\frac{\delta }{\beta }\frac{r_{2}S-p}{S+q}-\frac{mS+I}{k+N}-\beta I\\ +I-I_w^{\ast }\beta S-\delta \\ \le -\frac{r_1}{K}a+b{N}_w^{\ast }+c\frac{\delta }{\beta }{N-N_w^{\ast }}^2\\ +\frac{\theta a+b{N}_w^{\ast }}{a+bN+cS}\frac{{N-N_w^{\ast }}^2+{S-\delta /\beta }^2}{2}\\ +\frac{b\theta \delta }{\beta a+bN+cS}{N-N_w^{\ast }}^2\\ -\frac{m}{k+N}{S-\frac{\delta }{\beta }}^2+\frac{r_2{S}^2}{S+q}-\frac{p{r}_{2}S}{S+q}\\ -\frac{\delta }{\beta }\frac{r_{2}S}{S+q}+\frac{\delta }{\beta }\frac{p{r}_2}{S+q}-\frac{m\delta }{\beta k+N}S\\ +\frac{{\delta }^2}{{\beta }^2}\frac{m}{k+N}-\frac{m}{k+N}SI+\frac{\delta }{\beta }\frac{m}{k+N}I\\ -\beta S{I}_w^{\ast }+\delta I_w^{\ast }\\ \le -\frac{r_1}{K}a+b{N}_w^{\ast }-c\frac{\delta }{\beta }-\frac{\theta a+b{N}_w^{\ast }}{2a}\\ -\frac{b\theta \delta }{a\beta }{N-N_w^{\ast }}^2+\frac{\theta a+b{N}_w^{\ast }}{2a}\\ -\frac{m}{k+{\eta }_1+\epsilon }{S-\frac{\delta }{\beta }}^2+r_2\\ -\frac{p{r}_2}{q}-\beta I_w^{\ast }S+\frac{\delta p{r}_2}{\beta {\eta }_2+\epsilon +q}\\ +\frac{{\delta }^{2}m}{{\beta }^{2}k}+\frac{\delta m{\eta }_2+\epsilon }{\beta k}+\delta I_w^{\ast }.\end{array}$$

Since we have $\theta \le \text{min}{\mathrm{\Omega }}_1,{\mathrm{\Omega }}_2$ and ${\mathrm{\Omega }}_3\le -p\le {\mathrm{\Omega }}_4$, then for an $\epsilon >0$, we can choose the following condition:$$\begin{array}{c}\text{(79)}& \frac{r_1}{K}a+b{N}_w^{\ast }-c\frac{\delta }{\beta }\ge \frac{\theta a+b{N}_w^{\ast }}{2a}+\frac{b\theta \delta }{a\beta },\frac{\theta a+b{N}_w^{\ast }}{2a}\le \frac{m}{k+{\eta }_1+\epsilon },\\ r_2<\frac{p{r}_2}{q}+\beta I_w^{\ast },\text{\hspace{0.17em}and}\\ \frac{\delta p{r}_2}{\beta {\eta }_2+\epsilon +q}\le -\frac{{\delta }^{2}m}{{\beta }^{2}k}+\frac{\delta m{\eta }_2+\epsilon }{\beta k}\\ +\delta I_w^{\ast },\end{array}$$which implies ${}^C{{\displaystyle D}}_t^{\alpha }{L}_{5}N,S,I\le 0,\hspace{1em}\forall N,S,I\in {ℝ}_{+}^3$, and ${}^C{{\displaystyle D}}_t^{\alpha }{L}_{5}N,S,I=0$ at $E^{\ast }$. Based on that, the only invariant set on which ${}^C{{\displaystyle D}}_t^{\alpha }{L}_{5}N,S,I=0$ is the singleton $E^{\ast }$. Employing the same argument, it follows that $E^{\ast }$ is globally asymptotically stable.

4.2. For the Strong Allee Effect in Predator

Theorem 19.

If $p\ge \theta K+r_2{\eta }_3+q/r_2$, then ${\mathrm{\Xi }}_N=K,\mathrm{0,0}$ is globally asymptotically stable.

Proof.

Let $Nt,St,It$ be any positive solution of system (4). From the hypothesis of Theorem 12, an $\epsilon >0$ can be chosen such that $\theta K+r_2\le r_{2}p/{\eta }_3+\epsilon +q$. In the strong Allee effect in predator, we have $0\le St\le {\eta }_3+\epsilon$. Let us consider a positive definite Lyapunov function at ${\mathrm{\Xi }}_N=K,\mathrm{0,0}$ as follows:$$\begin{array}{}\text{(80)}& L_{6}N,S,I=N-K-K\ln \frac{N}{K}+S+I.\end{array}$$

By taking the fractional time derivative of $L_{6}N,S,I$ along the solutions of system (4) and using Lemma 6, one has the following equation:$$\begin{array}{c}\text{(81)}& {}^C{{\displaystyle D}}_t^{\alpha }{L}_{6}N,S,I\le \frac{N-K}{N}{}^C{{\displaystyle D}}_t^{\alpha }N+{}^C{{\displaystyle D}}_t^{\alpha }S+{}^C{{\displaystyle D}}_t^{\alpha }I=\frac{N-K}{N}{r}_{1}1-\frac{N}{K}-\frac{\theta S}{a+bN+cS}N\\ +S{r}_2\frac{S-p}{S+q}-\frac{mS+I}{k+N}-\delta I\\ \le -\frac{r_1}{K}{N-K}^2+\theta K+r_2-\frac{r_{2}p}{{\eta }_3+q}S\\ -\delta I.\end{array}$$

Regarding that $\theta K+r_2\le r_{2}p/{\eta }_3+\epsilon +q$, we have ${}^C{{\displaystyle D}}_t^{\alpha }{L}_{6}N,S,I\le 0$ for all $N,S,I\in {ℝ}_{+}^3$. Further, ${}^C{{\displaystyle D}}_t^{\alpha }{L}_{6}N,S,I=0$ implies that $N,S,I=K,\mathrm{0,0}$. Using Lemma 7, we conclude that $E^{\ast }$ is globally asymptotically stable.

5. Numerical Simulations

In this section, we demonstrate the numerical simulations based on the Adams–Bashforth–Moulton predictor-corrector method provided by Diethelm et al. [43] to verify the theoretical results established in the previous section. In addition, we also present the complex dynamics of system (4) such as the existence of bistability, forward, backward, saddle-node, and Hopf bifurcations as the effects of the disease and the fractional derivative.

The first simulation is given to observe the role of the transmission rate $\beta$ to the dynamics of system (4) for the weak Allee effect case. We consider a set of hypothetical parameters in (26), where some of the parameters are taken from [28].$$\begin{array}{c}\text{(82)}& r_1=0.5,K=5,\theta =0.65,a=1,b=1,c=0.1,r_2=0.1,k=1,m=0.01,p=-0.46,\\ \delta =0.2,q=3.\end{array}$$

In view of the existence conditions of the equilibrium points, Theorem 10 and Theorem 11, we plot a bifurcation diagram for the weak Allee effect case as shown in Figure 1. Here, $\alpha$ is set at $\alpha =0.952$ and varies the value of $\beta$ in the range of $\mathrm{0,0.4}$. With these parameter values, system (4) always has the equilibrium points $E_0$ and $E_N$ which is unstable according to Theorem 10 (a) and (b) and does not have a positive planar equilibrium point $E_{NS}$. Therefore, we only observe the dynamics of system (4) around $E_S$ (the solid and dash black curve), $E_{SI}$ (the solid and dash green curve), and $E^{\ast }$ (the solid and dash blue curve).

[figure(s) omitted; refer to PDF]

In Figure 1, when $\beta <{\beta }_1^{\ast }=0.0263$, the axial equilibrium point $E_S$ is the only local asymptotically stable equilibrium point of system (4). A phase portrait of this behavior is shown in Figure 2(a), for example, $\beta =0.01$. One can easily confirm that $0.5=r_1<\theta S_w/a+c{S}_w=2.8078$, $2.54=p+q<m/k{r}_2{S_w+q}^2=11.2463$, and $7.6049\approx S_w<\delta /\beta =20$, which satisfy the locally asymptotically stable condition for $E_S=\mathrm{0,7.60487},0$ as stated in Theorem 10(c). If we increase the value of parameter $\beta$ such that $\beta >{\beta }_1^{\ast }=0.0263$, $E_S$ loses its stability and the stable planar equilibrium point $E_{SI}$ arises via forward bifurcation. We confirm this behavior by numerical solutions in Figure 2(b). For $\beta =0.05$, we have the following conditions: $0.5=r_1<\delta \theta /a\beta +c\delta =1.8571$ and $0.10367=\beta r_{2}p+q/{\delta +\beta q}^2<m/\beta k=0.2<r_2\delta -\beta p/\delta \delta +\beta q=0.31857$ that fulfill the stability conditions of $E_{SI}=\mathrm{0,4,0.39524}$ in Theorem 10$di$ but do not satisfy the stability conditions of $E_S$ in Theorem 10$\text{c}$. Since our numerical simulation set $\alpha =0.952$, the stability of $E_{SI}$ is only attained at interval $\beta \in {\beta }_1^{\ast },{\beta }_4^{\ast }$, where ${\beta }_4^{\ast }\approx 0.24$. From the existence conditions of the interior point, we observe that when $0<\beta <{\beta }_2^{\ast }\approx 0.1377$, the interior equilibrium point $E^{\ast }$ does not exist. Next, based on the existence condition of the interior equilibrium point, two positive interior equilibrium points $E_1^{\ast }$ and $E_2^{\ast }$ which have different sign appear simultaneously after $\beta$ across the critical value ${\beta }_2^{\ast }\approx 0.1377$ via saddle-node bifurcation. Notice that for the set of parameters and ${\beta }_2^{\ast }<\beta <{\beta }_3^{\ast }$, the stability conditions of $E_{SI}$ and $E^{\ast }$ are satisfied. We plot a phase portrait for $\beta =0.14$ in Figure 2(c) to present the behavior, where $E_1^{\ast }=2.31322,1.42857,0.26803$ and $E_{SI}=\mathrm{0,1.428571429},0.1890629800$ are stable equilibrium points and $E_2^{\ast }=1.54392,1.42857,0.25727$ is the unstable equilibrium point. It means system (4) experiences the bistability phenomenon in that interval.

[figure(s) omitted; refer to PDF]

A further investigation shows that the stable branch interior point $E_1^{\ast }$ experiences a Hopf bifurcation when we increase $\beta$ passes through ${\beta }_3^{\ast }\approx 0.17475$. If we check the transversality condition through Figure 3(a), then we have ${d\chi \beta /d\beta }_{\beta ={\beta }_3^{\ast }}\approx 0.15602\pi \ne 0$. To confirm Hopf bifurcation, we choose $\beta =0.20$ (see Figure 2(d)). One can easily check that the interior point $E_1^{\ast }=3.62407,\mathrm{1,0.16985}$ has one real negative eigen value, a pair of complex conjugate eigen values with positive real part and $\alpha >{\alpha }^{\ast }=0.95044$. Therefore, $E_1^{\ast }$ fail to keep its stability which is indicated by the existence of limit cycle enclosing the interior point $E_1^{\ast }$. Similar behavior around the interior point $E_1^{\ast }$ still occur until $\beta$ reaches ${\beta }_5^{\ast }\approx 0.26407$ where the interior equilibrium point $E_1^{\ast }$ retrieves again its stability via Hopf bifurcation for $\beta$ larger than ${\beta }_5^{\ast }\approx 0.26407$. Here, we verify the transversality condition using Figure 3(b) which satisfies ${d\chi \beta /d\beta }_{\beta ={\beta }_5^{\ast }}\approx 0.01915\pi \ne 0$. Therefore, to illustrate this behavior, we depict a phase portrait of system (4) for $\beta =0.4$ in Figure 2(f). On the other hand, the unstable branch interior point $E_2^{\ast }$ collides with the stable planar equilibrium point $E_{SI}$ when $\beta$ attains ${\beta }_4^{\ast }$ via backward bifurcation. Hence, the planar equilibrium point $E_{SI}$ also loses its stability, as shown in Figure 2(e). In this case, all existing equilibrium points of system (4) lose their stability. However, as we can see in Figure 2(e), all solutions of system (4) converge to a limit cycle around the interior point $E_1^{\ast }$.

[figure(s) omitted; refer to PDF]

Next, we discuss the contribution of transmission rate $\beta$ to the dynamics of system (4) under the strong Alee effect case $p=0.4>0$ with the same set of parameter values as in [28] subsection 6.2.$$\begin{array}{c}\text{(83)}& r_1=0.5,K=5,\theta =0.4,a=1,b=1,c=0.1,r_2=0.1,k=1,m=0.05,p=0.4,\\ \delta =0.2,q=3,\alpha =\mathrm{0.98.}\end{array}$$

Notice that, for the strong Alee effect case $p>0$, ${\mathrm{\Xi }}_N=\mathrm{5,0,0}$ is always locally asymptotically stable according to Theorem 12$\text{b}$ and ${\mathrm{\Xi }}_0=\mathrm{0,0,0}$ is always a saddle point as stated in Theorem 12(a). In Figure 4(a), we portray a phase portrait of system (4) for $\beta =0.073$. System (4) has the following equilibrium points: ${\mathrm{\Xi }}_0=\mathrm{0,0,0},{\mathrm{\Xi }}_N=\mathrm{5,0,0},{\mathrm{\Xi }}_{N{S}_1}=4.58193,0.58956,0$, and ${\mathrm{\Xi }}_{N{S}_2}=2.2404,2.40121,0$. Since ${\mathrm{\Xi }}_{N{S}_1}$ satisfies the stability conditions in Theorem 12(e) (i), then system (4) has bistability phenomenon. When we increase $\beta$ to $\beta =0.10$, the planar equilibrium point ${\mathrm{\Xi }}_{N{S}_2}=2.2404,2.40121,0$ becomes a saddle point along with the emergence of the locally asymptotically stable interior point ${\mathrm{\Xi }}^{\ast }=3.16886,\mathrm{2,0.07154}$. In this case, we still have bistability condition, as we can see in Figure 4(b). Further, when we keep increasing $\beta$ until $\beta =0.73$, the interior point ${\mathrm{\Xi }}^{\ast }$ of system (4) does not exists anymore and the equilibrium point ${\mathrm{\Xi }}_N=\mathrm{5,0,0}$ becomes a unique stable equilibrium point of system (4), as shown in Figure 4(c).

[figure(s) omitted; refer to PDF]

Remark 20.

According to Figures 1, 2, and 4, the transmission rate of disease $\beta$ has a great impact to the dynamics of system (4). The numerical simulations indicate that system (4) could exhibit the forward, backward, saddle-node, and Hopf bifurcations driven by $\beta$. System (4) also performs a bi-stability phenomenon both for weak and strong Alee effects which suggest that solutions of system (4) are sensitive enough to the initial conditions. Figures 2(a), 2(b), and 2(f) confirm the global behavior of system (4) for the weak Allee effect.

Remark 21.

From an ecological point of view, as the transmission rate of disease $\beta$ increases in the predator population, the number of susceptible predator population decreases; on the contrary, the prey population increases. When the Allee effect is weak, there is no possibility of extinction for the predator, and conversely when the Allee effect is strong, there is always a possibility of extinction for the predator. We also found for the same Allee threshold, models with disease would have a higher risk of predator extinction when compared to models without disease. This result suggests that to save all of the species from extinction with disease threatening one of the population, one should control the transmission rate of disease.

In the next simulation, we will show the influence of the order of fractional derivative $\alpha$ on the dynamics of system (4) using the following parameter values:$$\begin{array}{c}\text{(84)}& r_1=0.5,K=5,\theta =0.4,a=1,b=1,c=0.1,r_2=0.1,k=1,m=0.05,p=-0.46,\\ \delta =0.2,\beta =0.073\\ q=3.\end{array}$$

For the above parameter values, it is found that system (4) has four equilibrium points, where $E_0=\mathrm{0,0,0},E_N=\mathrm{5,0,0},E_S=\mathrm{0,0.58166},0$ are saddle points, respectively, and $E_{NS}=1.21751,2.31593,0$ is a conditionally locally asymptotically stable equilibrium point. It can be observed that system (4) undergoes a Hopf bifurcation at ${\alpha }^{\ast }=0.84124$, as shown in Figure 5(a). We also can examine for $\alpha ={\alpha }^{\ast }$ the transversality condition is holds, i.e., ${\partial \chi \alpha /\partial \alpha |}_{\alpha ={\alpha }^{\ast }}=\pi /2\ne 0$. Then, we can observe that when $\alpha <{\alpha }^{\ast }=0.84124$, the planar equilibrium point $E_{NS}=1.21751,2.31593,0$ is locally asymptotically stable as exhibited in Figure 5(b) for $\alpha =0.82$. Otherwise, when $\alpha >{\alpha }^{\ast }$, there exists a limit cycle around $E_{NS}=1.21751,2.31593,0$ which shows that both prey and susceptible predator populations are fluctuating, as shown in Figure 5(c).

We next exhibit the existence of two limit cycles as the solutions of system (4) by varying the order of the derivative $\alpha$ at $\alpha =0.85,0.96$, and $0.98$. Using the set of parameters in (26), except parameter $\beta$ fixed at $\beta =0.22$, system (4) has six equilibrium points, i.e., $E_0=\mathrm{0,0,0}$, $E_N=\mathrm{5,0,0}$, $E_S=\mathrm{0,7.60487},0$, $E_{SI}=\mathrm{0,0.90909},0.11275$, $E_1^{\ast }=3.78913,0.90909,0.14915$, and $E_2^{\ast }=0.11996,0.90909,0.11753$. One can easily check that $E_0,E_N,E_S$, and $E_2^{\ast }$ are unstable, while the stability of $E_{SI}$ and $E_1^{\ast }$ depends on the order of the derivative $\alpha$. When $\alpha =0.85$, system (4) undergoes the bistability phenomenon which occurs around the interior equilibrium point $E_1^{\ast }$ and the planar equilibrium point $E_{SI}$ as we can see in Figure 6(a). This is confirmed by the real negative eigen value of $E_1^{\ast }$ and $E_{SI}$, also the critical value of the equilibrium points ${\alpha }^{\ast }$ are $0.9505$ and $0.9729$, respectively, which are greater than $\alpha =0.85$. The solutions of system (4) with two close enough initial conditions converge to different equilibrium points. When we increase $\alpha$ to $0.96$, $E_{SI}$ is still locally asymptotically stable $\text{since\hspace{0.17em}}0.9729={\alpha }^{\ast }>\alpha =0.96$, while $E_1^{\ast }$ becomes unstable via a Hopf bifurcation and limit cycle appears enclosing the interior point $E_1^{\ast }$, as shown in Figure 6(b). When we continue to increase $\alpha$ at $\alpha =0.98$, following $E_1^{\ast }$, then $E_{SI}$ loses its stability too through a Hopf bifurcation and two stable limit cycles appear instead in system (4), see Figure 6(c). These numerical simulations show that system (4) could exhibit two Hopf bifurcations simultaneously which is controlled by the order of the derivative $\alpha$.

[figure(s) omitted; refer to PDF]

Remark 22.

From the ecological point of view, a phenomenon which is given in Figure 6(c) show that if the initial value of prey population is relatively small then the susceptible and infected predators will oscillate even in the absence of prey population. On the other side, when the initial prey population is relatively large, all populations may oscillate for a long period of time.

Remark 23.

The analytical results and numerical simulations of the model in the absence of disease in predator population (system (1)) had been found by Rahmi et al. [28]. Numerically, the recent model with disease in predator population has more rich dynamical behaviors than the former model. One of the reasons comes from varying the transmission rate of disease parameter. When there is no disease, there is only bistability, forward and Hopf bifurcations, meanwhile there is bistability, saddle-node, backward, forward, and Hopf bifurcations for the model with disease. More specifically, the proposed model exhibits a new phenomenon in the fractional order system, which is the existence of two limit cycles driven by the order of derivative as given in Figure 6(c).

6. Conclusion

In this paper, we have merged a predator-prey model and epidemiology model into the eco-epidemiological model. In ecological literature, there is much evidence that the double Allee effect may be acting on a single population. Moreover, it is natural to assume that there is a deadly infectious disease spreading within the population. Thus, we have developed an eco-epidemiological model incorporating the double Allee effect and disease spread on the predator population. This eco-epidemiological model has been modelled by a system with fractional-order differential equations in Caputo sense. We have split the predator population into two subpopulations: the susceptible predator and the infected predator. We proved that the solution of system (4) exists uniquely and whenever we started with a positive initial condition then all solutions remain non-negative and bounded. We showed that our system (4) has four types of equilibrium points, i.e., the trivial, two axial, two planar, and the interior equilibrium points. The trivial equilibrium point for both the strong and weak Allee effect cases is always unstable, which means that there is no condition of extinction of all populations in the future. The first axial equilibrium point (the predator extinction point) is always unstable for the weak Allee effect case, while for the strong Alee effect case, it is always stable. It means that the predator population can go to extinction in the strong Allee effect, which is contrary to the weak Allee effect case. The second axial equilibrium point (the prey and infected predator extinction point) is conditionally stable. We have observed that the dynamics of the axial equilibrium point do not depend on the order of the fractional derivative. Then, we have two planar equilibrium points (the prey extinction point and the infected predator extinction point) which are conditionally stable. The interior equilibrium point is also conditionally stable for both weak and strong Allee effects. The stability conditions of the last two types of equilibrium points show that the order of the fractional derivative $\alpha$ could affect its stability. By choosing an order $\alpha$ which is smaller than the critical value of ${\alpha }^{\ast }$, then the equilibrium point may be stabilized. Numerical simulations show that the dynamics of system (4) have corresponded well with our theoretical results. It has been observed numerically that the transmission rate $\beta$ takes a vital role to control the dynamical behavior of system (4). The combination of disease and double Allee effects on predators may accelerate the decline in predator populations which ultimately increases the risk of extinction. We also show numerically the existence of the Hopf, forward, backward, and saddle-node bifurcations driven by the transmission rate $\beta$. Further, the bistability phenomenon also appears in both the strong and weak Allee effect cases. All numerical results indicate that each population in system (4) could be extinct or survive or oscillate depending on the parameter values and initial conditions. In [28], the authors studied the dynamics of the model with no disease in the predator. They also carried out some numerical simulations to show how the capturing rate, Allee threshold, and the order of the derivative affect the dynamics of the model. The comparative results of these two models, with and without disease in the predator population are given in Remark 23. However, it could be interesting if the eco-epidemiological model with the double Allee effect and disease in predators also considered vaccination strategies as the approach to reduce population extinction risks. For example, there are African wild dogs that are endangered carnivores which get vaccination coverage for rabies virus and canine distemper virus (CDV) [44].

Authors’ Contributions

All authors have contributed equally and have finalized the manuscript.

Acknowledgments

This work was supported by the FMIPA via DPA PTNBH-University of Brawijaya, according to the Doctoral (Associate Professor) Research Grant, with contract number: 3110.12/UN10.F09/PN/2022, dated May 30, 2022.

References

[1] P. J. Pal, T. Saha, "Qualitative analysis of a predator-prey system with double Allee effect in prey," Chaos, Solitons & Fractals, vol. 73, pp. 36-63, DOI: 10.1016/j.chaos.2014.12.007, 2015.

[2] Y. Vera-Damián, C. Vidal, E. González-Olivares, "Dynamics and bifurcations of a modified Leslie-Gower–type model considering a Beddington-DeAngelis functional response," Mathematical Methods in the Applied Sciences, vol. 42 no. 9, pp. 3179-3210, 2019.

[3] S. Yu, "Global stability of a modified Leslie-Gower model with Beddington-DeAngelis functional response," Advances in Difference Equations, vol. 2014 no. 1,DOI: 10.1186/1687-1847-2014-84, 2014.

[4] S. Vinoth, R. Sivasamy, K. Sathiyanathan, B. Unyong, G. Rajchakit, R. Vadivel, N. Gunasekaran, "The dynamics of a Leslie type predator-prey model with fear and Allee effect, " , Advances in Different Equations, vol. 798 no. 1, 2001.

[5] S. Vinoth, R. Sivasamy, K. Sathiyanathan, B. Unyong, G. Vadivel, R. Gunasekaran, N. Gunasekaran, "A novel discrete-time leslie-gower model with the impact of allee effect in predator population," Complexity, vol. 2022,DOI: 10.1155/2022/6931354, 2022.

[6] A. Suryanto, I. Darti, "Dynamics of Leslie-Gower pest-predator model with disease in pest including pest-harvesting and optimal implementation of pesticide," International Journal of Mathematics and Mathematical Sciences, vol. 2019,DOI: 10.1155/2019/5079171, 2019.

[7] S. Sharma, G. P. Samanta, "A Leslie-Gower predator-prey model with disease in prey incorporating a prey refuge," Chaos, Solitons & Fractals, vol. 70 no. 1, pp. 69-84, DOI: 10.1016/j.chaos.2014.11.010, 2015.

[8] S. Wang, Z. Ma, W. Wang, "Dynamical behavior of a generalized eco-epidemiological system with prey refuge," Advances in Difference Equations, vol. 2018 no. 1,DOI: 10.1186/s13662-018-1704-x, 2018.

[9] B. Ghanbari, S. Djilali, "Mathematical and numerical analysis of a three‐species predator‐prey model with herd behavior and time fractional‐order derivative," Mathematical Methods in the Applied Sciences, vol. 43 no. 4, pp. 1736-1752, DOI: 10.1002/mma.5999, 2020.

[10] U. Kumar, P. S. Mandal, E. Venturino, "Impact of Allee effect on an eco-epidemiological system," Ecological Complexity, vol. 42 no. 100828,DOI: 10.1016/j.ecocom.2020.100828, 2020.

[11] W. C. Allee, Animal Aggregations, a Study in General Sociology, 1931.

[12] F. Courchamp, L. Berec, J. Gascoigne, Allee Effects in Ecology and Conservation, 2008.

[13] E. Rahmi, I. Darti, A. Suryanto, Trisilowati, H. S. Panigoro, "Stability analysis of a fractional-order leslie-gower model with allee effect in predator," Journal of Physics: Conference Series, vol. 1821 no. 012051,DOI: 10.1088/1742-6596/1821/1/012051, 2021.

[14] S. K. Sasmal, D. S. Mandal, J. Chattopadhyay, "A predator-pest model with Allee effect and pest culling and additional food provision to the predator - application to pest control," Journal of Biological Systems, vol. 25 no. 02, pp. 295-326, DOI: 10.1142/s0218339017500152, 2017.

[15] A. Rojas-Palma, E. González-Olivares, "Optimal harvesting in a predator-prey model with Allee effect and sigmoid functional response," Applied Mathematical Modelling, vol. 36 no. 5, pp. 1864-1874, DOI: 10.1016/j.apm.2011.07.081, 2012.

[16] P. D. Srinivasu, G. Kiran Kumar, "Bio-economics of a renewable resource subjected to strong Allee effect," Communications in Nonlinear Science and Numerical Simulation, vol. 19 no. 6, pp. 1686-1696, DOI: 10.1016/j.cnsns.2013.10.016, 2014.

[17] M. K. Singh, B. S. Bhadauria, B. K. Singh, "Bifurcation analysis of modified Leslie-Gower predator-prey model with double Allee effect," Ain Shams Engineering Journal, vol. 9 no. 4, pp. 1263-1277, DOI: 10.1016/j.asej.2016.07.007, 2018.

[18] L. Berec, E. Angulo, F. Courchamp, "Multiple Allee effects and population management," Trends in Ecology & Evolution, vol. 22 no. 4, pp. 185-191, DOI: 10.1016/j.tree.2006.12.002, 2007.

[19] M. S. Abdelouahab, N. E. Hamri, J. Wang, "Hopf bifurcation and chaos in fractional-order modified hybrid optical system," Nonlinear Dynamics, vol. 69 no. 1–2, pp. 275-284, DOI: 10.1007/s11071-011-0263-4, 2012.

[20] P. Huang, J. Zhang, N. Jean Damascene, C. Dong, Z. Wang, "A fractional order viscoelastic-plastic creep model for coal sample considering initial damage accumulation," Alexandria Engineering Journal, vol. 60 no. 4, pp. 3921-3930, DOI: 10.1016/j.aej.2021.02.054, 2021.

[21] S. Daudi, L. Luboobi, M. Kgosimore, D. Kuznetsov, "A fractional-order fall armyworm-maize biomass model with naturally beneficial insects and optimal farming awareness," Results in Applied Mathematics, vol. 12 no. 2021,DOI: 10.1016/j.rinam.2021.100209, 2021.

[22] T. B. Gashirai, S. D. Musekwa-Hove, P. O. Lolika, S. Mushayabasa, "Global stability and optimal control analysis of a foot-and-mouth disease model with vaccine failure and environmental transmission," Chaos, Solitons & Fractals, vol. 132 no. 2020,DOI: 10.1016/j.chaos.2019.109568, 2020.

[23] J. Alidousti, E. Ghafari, "Dynamic behavior of a fractional order prey-predator model with group defense," Chaos, Solitons & Fractals, vol. 134 no. 2020,DOI: 10.1016/j.chaos.2020.109688, 2020.

[24] M. Das, G. P. Samanta, "A prey-predator fractional order model with fear effect and group defense," International Journal of Dynamics and Control, vol. 9 no. 1, pp. 334-349, DOI: 10.1007/s40435-020-00626-x, 2021.

[25] H. S. Panigoro, A. Suryanto, W. M. Kusumawinahyu, I. Darti, "A Rosenzweig–MacArthur model with continuous threshold harvesting in predator involving fractional derivatives with power law and Mittag–Leffler kernel," Axioms, vol. 9 no. 4,DOI: 10.3390/axioms9040122, 2020.

[26] K. Diethelm, "A fractional calculus based model for the simulation of an outbreak of dengue fever," Nonlinear Dynamics, vol. 71 no. 4, pp. 613-619, DOI: 10.1007/s11071-012-0475-2, 2013.

[27] R. Almeida, N. R. O. Bastos, M. T. T. Monteiro, "Modeling some real phenomena by fractional differential equations," Mathematical Methods in the Applied Sciences, vol. 39 no. 16, pp. 4846-4855, DOI: 10.1002/mma.3818, 2016.

[28] E. Rahmi, I. Darti, A. Suryanto, Trisilowati, "A modified leslie-gower model incorporating beddington-DeAngelis functional response, double allee effect and memory effect," Fractal and Fractional, vol. 5 no. 3,DOI: 10.3390/fractalfract5030084, 2021.

[29] J. R. Beddington, "Mutual interference between parasites or predators and its effect on searching efficiency," Journal of Animal Ecology, vol. 44 no. 1,DOI: 10.2307/3866, 1975.

[30] A. Mondal, A. K. Pal, G. P. Samanta, "On the dynamics of evolutionary Leslie-Gower predator-prey eco-epidemiological model with disease in predator," Ecological Genetics and Genomics, vol. 10,DOI: 10.1016/j.egg.2018.11.002, 2019.

[31] A. Kang, Y. Xue, J. Fu, "Dynamic behaviors of a Leslie-Gower ecoepidemiological model," Discrete Dynamics in Nature and Society, vol. 2015,DOI: 10.1155/2015/169242, 2015.

[32] I. Petras, Fractional-Order Nonlinear Systems, 2011.

[33] Z. M. Odibat, N. T. Shawagfeh, "Generalized Taylor's formula," Applied Mathematics and Computation, vol. 186 no. 1, pp. 286-293, DOI: 10.1016/j.amc.2006.07.102, 2007.

[34] H.-L. Li, L. Zhang, C. Hu, Y.-L. Jiang, Z. Teng, "Dynamical analysis of a fractional-order predator-prey model incorporating a prey refuge," Journal of Applied Mathematics and Computing, vol. 54 no. 1–2, pp. 435-449, DOI: 10.1007/s12190-016-1017-8, 2017.

[35] D. Matignon, "Stability results for fractional differential equations with applications to control processing," Computational Engineering in Systems Applications, pp. 963-968, 1996.

[36] C. Vargas-De-León, "Volterra-type Lyapunov functions for fractional-order epidemic systems," Communications in Nonlinear Science and Numerical Simulation, vol. 24 no. 1–3, pp. 75-85, DOI: 10.1016/j.cnsns.2014.12.013, 2015.

[37] J. Huo, H. Zhao, L. Zhu, "The effect of vaccines on backward bifurcation in a fractional order HIV model," Nonlinear Analysis: Real World Applications, vol. 26, pp. 289-305, DOI: 10.1016/j.nonrwa.2015.05.014, 2015.

[38] S. K. Choi, B. Kang, N. Koo, "Stability for Caputo fractional differential systems," Abstract and Applied Analysis, vol. 2014,DOI: 10.1155/2014/631419, 2014.

[39] E. Ahmed, A. El-Sayed, H. A. El-Saka, "On some Routh-Hurwitz conditions for fractional order differential equations and their applications in Lorenz, Rössler, Chua and Chen systems," Physics Letters A, vol. 358 no. 1,DOI: 10.1016/j.physleta.2006.04.087, 2006.

[40] H. S. Panigoro, A. Suryanto, W. M. Kusumawinahyu, I. Darti, "Dynamics of an eco-epidemic predator-prey model involving fractional derivatives with power-law and mittag-leffler kernel," Symmetry, vol. 13 no. 5,DOI: 10.3390/sym13050785, 2021.

[41] M. S. Tavazoei, "A note on fractional-order derivatives of periodic functions," Automatica, vol. 46 no. 5, pp. 945-948, DOI: 10.1016/j.automatica.2010.02.023, 2010.

[42] X. Li, R. Wu, "Hopf bifurcation analysis of a new commensurate fractional-order hyperchaotic system," Nonlinear Dynamics, vol. 78 no. 1, pp. 279-288, DOI: 10.1007/s11071-014-1439-5, 2014.

[43] K. Diethelm, N. J. Ford, A. D. Freed, "A predictor-corrector approach for the numerical solution of fractional differential equations," Nonlinear Dynamics, vol. 29 no. 1–4,DOI: 10.1023/a:1016592219341, 2002.

[44] K. C. Prager, R. Woodroffe, A. Cameron, D. T. Haydon, "Vaccination strategies to conserve the endangered African wild dog (Lycaon pictus)," Biological Conservation, vol. 144 no. 7, pp. 1940-1948, DOI: 10.1016/j.biocon.2011.03.025, 2011.